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