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 (at) 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 94
