1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Gael Guennebaud <[email protected]> 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 EIGEN_DEVICE_FUNC 28 void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha) 29 { 30 const Index size = u.size(); 31 for (Index i=0; i<size; ++i) 32 { 33 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) += 34 (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.tail(size-i) 35 + (alpha * numext::conj(v.coeff(i))) * u.tail(size-i); 36 } 37 } 38 }; 39 40 template<typename Scalar, typename Index, typename UType, typename VType> 41 struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Upper> 42 { 43 static void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha) 44 { 45 const Index size = u.size(); 46 for (Index i=0; i<size; ++i) 47 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) += 48 (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.head(i+1) 49 + (alpha * numext::conj(v.coeff(i))) * u.head(i+1); 50 } 51 }; 52 53 template<bool Cond, typename T> struct conj_expr_if 54 : conditional<!Cond, const T&, 55 CwiseUnaryOp<scalar_conjugate_op<typename traits<T>::Scalar>,T> > {}; 56 57 } // end namespace internal 58 59 template<typename MatrixType, unsigned int UpLo> 60 template<typename DerivedU, typename DerivedV> 61 EIGEN_DEVICE_FUNC SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> 62 ::rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha) 63 { 64 typedef internal::blas_traits<DerivedU> UBlasTraits; 65 typedef typename UBlasTraits::DirectLinearAccessType ActualUType; 66 typedef typename internal::remove_all<ActualUType>::type _ActualUType; 67 typename internal::add_const_on_value_type<ActualUType>::type actualU = UBlasTraits::extract(u.derived()); 68 69 typedef internal::blas_traits<DerivedV> VBlasTraits; 70 typedef typename VBlasTraits::DirectLinearAccessType ActualVType; 71 typedef typename internal::remove_all<ActualVType>::type _ActualVType; 72 typename internal::add_const_on_value_type<ActualVType>::type actualV = VBlasTraits::extract(v.derived()); 73 74 // If MatrixType is row major, then we use the routine for lower triangular in the upper triangular case and 75 // vice versa, and take the complex conjugate of all coefficients and vector entries. 76 77 enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 }; 78 Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived()) 79 * numext::conj(VBlasTraits::extractScalarFactor(v.derived())); 80 if (IsRowMajor) 81 actualAlpha = numext::conj(actualAlpha); 82 83 typedef typename internal::remove_all<typename internal::conj_expr_if<int(IsRowMajor) ^ int(UBlasTraits::NeedToConjugate), _ActualUType>::type>::type UType; 84 typedef typename internal::remove_all<typename internal::conj_expr_if<int(IsRowMajor) ^ int(VBlasTraits::NeedToConjugate), _ActualVType>::type>::type VType; 85 internal::selfadjoint_rank2_update_selector<Scalar, Index, UType, VType, 86 (IsRowMajor ? int(UpLo==Upper ? Lower : Upper) : UpLo)> 87 ::run(_expression().const_cast_derived().data(),_expression().outerStride(),UType(actualU),VType(actualV),actualAlpha); 88 89 return *this; 90 } 91 92 } // end namespace Eigen 93 94 #endif // EIGEN_SELFADJOINTRANK2UPTADE_H 95