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, 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 | (conj(alpha) * conj(u.coeff(i))) * v.tail(size-i) |
---|
34 | + (alpha * 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, 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 | (conj(alpha) * conj(u.coeff(i))) * v.head(i+1) |
---|
48 | + (alpha * 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, 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 | * internal::conj(VBlasTraits::extractScalarFactor(v.derived())); |
---|
79 | if (IsRowMajor) |
---|
80 | actualAlpha = internal::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 |
---|