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 13namespace Eigen { 14 15namespace internal { 16 17/* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu' 18 * It corresponds to the Level2 syr2 BLAS routine 19 */ 20 21template<typename Scalar, typename Index, typename UType, typename VType, int UpLo> 22struct selfadjoint_rank2_update_selector; 23 24template<typename Scalar, typename Index, typename UType, typename VType> 25struct 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 39template<typename Scalar, typename Index, typename UType, typename VType> 40struct 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 52template<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 58template<typename MatrixType, unsigned int UpLo> 59template<typename DerivedU, typename DerivedV> 60SelfAdjointView<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