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 |
---|