[136] | 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
|
---|