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_SELFADJOINT_PRODUCT_H 11 #define EIGEN_SELFADJOINT_PRODUCT_H 12 13 /********************************************************************** 14 * This file implements a self adjoint product: C += A A^T updating only 15 * half of the selfadjoint matrix C. 16 * It corresponds to the level 3 SYRK and level 2 SYR Blas routines. 17 **********************************************************************/ 18 19 namespace Eigen { 20 21 22 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> 23 struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs> 24 { 25 static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha) 26 { 27 internal::conj_if<ConjRhs> cj; 28 typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; 29 typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjLhsType; 30 for (Index i=0; i<size; ++i) 31 { 32 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1))) 33 += (alpha * cj(vecY[i])) * ConjLhsType(OtherMap(vecX+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1))); 34 } 35 } 36 }; 37 38 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> 39 struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs> 40 { 41 static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha) 42 { 43 selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vecY,vecX,alpha); 44 } 45 }; 46 47 template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime> 48 struct selfadjoint_product_selector; 49 50 template<typename MatrixType, typename OtherType, int UpLo> 51 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true> 52 { 53 static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha) 54 { 55 typedef typename MatrixType::Scalar Scalar; 56 typedef typename MatrixType::Index Index; 57 typedef internal::blas_traits<OtherType> OtherBlasTraits; 58 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType; 59 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType; 60 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived()); 61 62 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived()); 63 64 enum { 65 StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor, 66 UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1 67 }; 68 internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other; 69 70 ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(), 71 (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data())); 72 73 if(!UseOtherDirectly) 74 Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther; 75 76 selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo, 77 OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex, 78 (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex> 79 ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualOtherPtr, actualAlpha); 80 } 81 }; 82 83 template<typename MatrixType, typename OtherType, int UpLo> 84 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false> 85 { 86 static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha) 87 { 88 typedef typename MatrixType::Scalar Scalar; 89 typedef typename MatrixType::Index Index; 90 typedef internal::blas_traits<OtherType> OtherBlasTraits; 91 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType; 92 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType; 93 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived()); 94 95 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived()); 96 97 enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 }; 98 99 internal::general_matrix_matrix_triangular_product<Index, 100 Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex, 101 Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex, 102 MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo> 103 ::run(mat.cols(), actualOther.cols(), 104 &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(), 105 mat.data(), mat.outerStride(), actualAlpha); 106 } 107 }; 108 109 // high level API 110 111 template<typename MatrixType, unsigned int UpLo> 112 template<typename DerivedU> 113 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> 114 ::rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha) 115 { 116 selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha); 117 118 return *this; 119 } 120 121 } // end namespace Eigen 122 123 #endif // EIGEN_SELFADJOINT_PRODUCT_H 124
