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_SELFADJOINTRANK2UPTADE_H
|
---|
11 | #define EIGEN_SELFADJOINTRANK2UPTADE_H
|
---|
12 |
|
---|
13 | namespace Eigen {
|
---|
14 |
|
---|
15 | namespace internal {
|
---|
16 |
|
---|
17 | /* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu'
|
---|
18 | * It corresponds to the Level2 syr2 BLAS routine
|
---|
19 | */
|
---|
20 |
|
---|
21 | template<typename Scalar, typename Index, typename UType, typename VType, int UpLo>
|
---|
22 | struct selfadjoint_rank2_update_selector;
|
---|
23 |
|
---|
24 | template<typename Scalar, typename Index, typename UType, typename VType>
|
---|
25 | struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Lower>
|
---|
26 | {
|
---|
27 | static void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha)
|
---|
28 | {
|
---|
29 | const Index size = u.size();
|
---|
30 | for (Index i=0; i<size; ++i)
|
---|
31 | {
|
---|
32 | Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) +=
|
---|
33 | (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.tail(size-i)
|
---|
34 | + (alpha * numext::conj(v.coeff(i))) * u.tail(size-i);
|
---|
35 | }
|
---|
36 | }
|
---|
37 | };
|
---|
38 |
|
---|
39 | template<typename Scalar, typename Index, typename UType, typename VType>
|
---|
40 | struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Upper>
|
---|
41 | {
|
---|
42 | static void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha)
|
---|
43 | {
|
---|
44 | const Index size = u.size();
|
---|
45 | for (Index i=0; i<size; ++i)
|
---|
46 | Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) +=
|
---|
47 | (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.head(i+1)
|
---|
48 | + (alpha * numext::conj(v.coeff(i))) * u.head(i+1);
|
---|
49 | }
|
---|
50 | };
|
---|
51 |
|
---|
52 | template<bool Cond, typename T> struct conj_expr_if
|
---|
53 | : conditional<!Cond, const T&,
|
---|
54 | CwiseUnaryOp<scalar_conjugate_op<typename traits<T>::Scalar>,T> > {};
|
---|
55 |
|
---|
56 | } // end namespace internal
|
---|
57 |
|
---|
58 | template<typename MatrixType, unsigned int UpLo>
|
---|
59 | template<typename DerivedU, typename DerivedV>
|
---|
60 | SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
|
---|
61 | ::rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha)
|
---|
62 | {
|
---|
63 | typedef internal::blas_traits<DerivedU> UBlasTraits;
|
---|
64 | typedef typename UBlasTraits::DirectLinearAccessType ActualUType;
|
---|
65 | typedef typename internal::remove_all<ActualUType>::type _ActualUType;
|
---|
66 | typename internal::add_const_on_value_type<ActualUType>::type actualU = UBlasTraits::extract(u.derived());
|
---|
67 |
|
---|
68 | typedef internal::blas_traits<DerivedV> VBlasTraits;
|
---|
69 | typedef typename VBlasTraits::DirectLinearAccessType ActualVType;
|
---|
70 | typedef typename internal::remove_all<ActualVType>::type _ActualVType;
|
---|
71 | typename internal::add_const_on_value_type<ActualVType>::type actualV = VBlasTraits::extract(v.derived());
|
---|
72 |
|
---|
73 | // If MatrixType is row major, then we use the routine for lower triangular in the upper triangular case and
|
---|
74 | // vice versa, and take the complex conjugate of all coefficients and vector entries.
|
---|
75 |
|
---|
76 | enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
|
---|
77 | Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived())
|
---|
78 | * numext::conj(VBlasTraits::extractScalarFactor(v.derived()));
|
---|
79 | if (IsRowMajor)
|
---|
80 | actualAlpha = numext::conj(actualAlpha);
|
---|
81 |
|
---|
82 | internal::selfadjoint_rank2_update_selector<Scalar, Index,
|
---|
83 | typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::type>::type,
|
---|
84 | typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::type>::type,
|
---|
85 | (IsRowMajor ? int(UpLo==Upper ? Lower : Upper) : UpLo)>
|
---|
86 | ::run(_expression().const_cast_derived().data(),_expression().outerStride(),actualU,actualV,actualAlpha);
|
---|
87 |
|
---|
88 | return *this;
|
---|
89 | }
|
---|
90 |
|
---|
91 | } // end namespace Eigen
|
---|
92 |
|
---|
93 | #endif // EIGEN_SELFADJOINTRANK2UPTADE_H
|
---|