SelfadjointRank2Update.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_SELFADJOINTRANK2UPTADE_H
26 #define EIGEN_SELFADJOINTRANK2UPTADE_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 
32 /* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu'
33  * It corresponds to the Level2 syr2 BLAS routine
34  */
35 
36 template<typename Scalar, typename Index, typename UType, typename VType, int UpLo>
37 struct selfadjoint_rank2_update_selector;
38 
39 template<typename Scalar, typename Index, typename UType, typename VType>
40 struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Lower>
41 {
42  static void run(Scalar* mat, Index stride, const UType& u, const VType& v, Scalar alpha)
43  {
44  const Index size = u.size();
45  for (Index i=0; i<size; ++i)
46  {
47  Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) +=
48  (conj(alpha) * conj(u.coeff(i))) * v.tail(size-i)
49  + (alpha * conj(v.coeff(i))) * u.tail(size-i);
50  }
51  }
52 };
53 
54 template<typename Scalar, typename Index, typename UType, typename VType>
55 struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Upper>
56 {
57  static void run(Scalar* mat, Index stride, const UType& u, const VType& v, Scalar alpha)
58  {
59  const Index size = u.size();
60  for (Index i=0; i<size; ++i)
61  Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) +=
62  (conj(alpha) * conj(u.coeff(i))) * v.head(i+1)
63  + (alpha * conj(v.coeff(i))) * u.head(i+1);
64  }
65 };
66 
67 template<bool Cond, typename T> struct conj_expr_if
68  : conditional<!Cond, const T&,
69  CwiseUnaryOp<scalar_conjugate_op<typename traits<T>::Scalar>,T> > {};
70 
71 } // end namespace internal
72 
73 template<typename MatrixType, unsigned int UpLo>
74 template<typename DerivedU, typename DerivedV>
75 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
77 {
78  typedef internal::blas_traits<DerivedU> UBlasTraits;
79  typedef typename UBlasTraits::DirectLinearAccessType ActualUType;
80  typedef typename internal::remove_all<ActualUType>::type _ActualUType;
81  typename internal::add_const_on_value_type<ActualUType>::type actualU = UBlasTraits::extract(u.derived());
82 
83  typedef internal::blas_traits<DerivedV> VBlasTraits;
84  typedef typename VBlasTraits::DirectLinearAccessType ActualVType;
85  typedef typename internal::remove_all<ActualVType>::type _ActualVType;
86  typename internal::add_const_on_value_type<ActualVType>::type actualV = VBlasTraits::extract(v.derived());
87 
88  // If MatrixType is row major, then we use the routine for lower triangular in the upper triangular case and
89  // vice versa, and take the complex conjugate of all coefficients and vector entries.
90 
91  enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
92  Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived())
93  * internal::conj(VBlasTraits::extractScalarFactor(v.derived()));
94  if (IsRowMajor)
95  actualAlpha = internal::conj(actualAlpha);
96 
97  internal::selfadjoint_rank2_update_selector<Scalar, Index,
98  typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::type>::type,
99  typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::type>::type,
100  (IsRowMajor ? int(UpLo==Upper ? Lower : Upper) : UpLo)>
101  ::run(_expression().const_cast_derived().data(),_expression().outerStride(),actualU,actualV,actualAlpha);
102 
103  return *this;
104 }
105 
106 } // end namespace Eigen
107 
108 #endif // EIGEN_SELFADJOINTRANK2UPTADE_H