• Home
  • History
  • Annotate
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
22 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
23 struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs>
24 {
25   static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha)
26   {
27     internal::conj_if<ConjRhs> cj;
28     typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
29     typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjLhsType;
30     for (Index i=0; i<size; ++i)
31     {
32       Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1)))
33           += (alpha * cj(vecY[i])) * ConjLhsType(OtherMap(vecX+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1)));
34     }
35   }
36 };
37 
38 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
39 struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs>
40 {
41   static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha)
42   {
43     selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vecY,vecX,alpha);
44   }
45 };
46 
47 template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime>
48 struct selfadjoint_product_selector;
49 
50 template<typename MatrixType, typename OtherType, int UpLo>
51 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
52 {
53   static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha)
54   {
55     typedef typename MatrixType::Scalar Scalar;
56     typedef typename MatrixType::Index Index;
57     typedef internal::blas_traits<OtherType> OtherBlasTraits;
58     typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
59     typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
60     typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
61 
62     Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
63 
64     enum {
65       StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor,
66       UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
67     };
68     internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
69 
70     ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
71       (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
72 
73     if(!UseOtherDirectly)
74       Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
75 
76     selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
77                               OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
78                             (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
79           ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualOtherPtr, actualAlpha);
80   }
81 };
82 
83 template<typename MatrixType, typename OtherType, int UpLo>
84 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
85 {
86   static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha)
87   {
88     typedef typename MatrixType::Scalar Scalar;
89     typedef typename MatrixType::Index Index;
90     typedef internal::blas_traits<OtherType> OtherBlasTraits;
91     typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
92     typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
93     typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
94 
95     Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
96 
97     enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
98 
99     internal::general_matrix_matrix_triangular_product<Index,
100       Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor,   OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
101       Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
102       MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo>
103       ::run(mat.cols(), actualOther.cols(),
104             &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
105             mat.data(), mat.outerStride(), actualAlpha);
106   }
107 };
108 
109 // high level API
110 
111 template<typename MatrixType, unsigned int UpLo>
112 template<typename DerivedU>
113 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
114 ::rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha)
115 {
116   selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
117 
118   return *this;
119 }
120 
121 } // end namespace Eigen
122 
123 #endif // EIGEN_SELFADJOINT_PRODUCT_H
124