SelfadjointProduct.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_SELFADJOINT_PRODUCT_H
26 #define EIGEN_SELFADJOINT_PRODUCT_H
27 
28 /**********************************************************************
29 * This file implements a self adjoint product: C += A A^T updating only
30 * half of the selfadjoint matrix C.
31 * It corresponds to the level 3 SYRK and level 2 SYR Blas routines.
32 **********************************************************************/
33 
34 namespace Eigen {
35 
36 template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs>
37 struct selfadjoint_rank1_update;
38 
39 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
40 struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs>
41 {
42  static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha)
43  {
44  internal::conj_if<ConjRhs> cj;
45  typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
46  typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjRhsType;
47  for (Index i=0; i<size; ++i)
48  {
49  Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1)))
50  += (alpha * cj(vec[i])) * ConjRhsType(OtherMap(vec+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1)));
51  }
52  }
53 };
54 
55 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
56 struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs>
57 {
58  static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha)
59  {
60  selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vec,alpha);
61  }
62 };
63 
64 template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime>
65 struct selfadjoint_product_selector;
66 
67 template<typename MatrixType, typename OtherType, int UpLo>
68 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
69 {
70  static void run(MatrixType& mat, const OtherType& other, typename MatrixType::Scalar alpha)
71  {
72  typedef typename MatrixType::Scalar Scalar;
73  typedef typename MatrixType::Index Index;
74  typedef internal::blas_traits<OtherType> OtherBlasTraits;
75  typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
76  typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
77  typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
78 
79  Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
80 
81  enum {
82  StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor,
83  UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
84  };
85  internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
86 
87  ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
88  (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
89 
90  if(!UseOtherDirectly)
91  Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
92 
93  selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
94  OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
95  (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
96  ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha);
97  }
98 };
99 
100 template<typename MatrixType, typename OtherType, int UpLo>
101 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
102 {
103  static void run(MatrixType& mat, const OtherType& other, typename MatrixType::Scalar alpha)
104  {
105  typedef typename MatrixType::Scalar Scalar;
106  typedef typename MatrixType::Index Index;
107  typedef internal::blas_traits<OtherType> OtherBlasTraits;
108  typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
109  typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
110  typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
111 
112  Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
113 
114  enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
115 
116  internal::general_matrix_matrix_triangular_product<Index,
117  Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex,
118  Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
119  MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo>
120  ::run(mat.cols(), actualOther.cols(),
121  &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
122  mat.data(), mat.outerStride(), actualAlpha);
123  }
124 };
125 
126 // high level API
127 
128 template<typename MatrixType, unsigned int UpLo>
129 template<typename DerivedU>
132 {
133  selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
134 
135  return *this;
136 }
137 
138 } // end namespace Eigen
139 
140 #endif // EIGEN_SELFADJOINT_PRODUCT_H