1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_SELFADJOINT_PRODUCT_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_SELFADJOINT_PRODUCT_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/********************************************************************** 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* This file implements a self adjoint product: C += A A^T updating only 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* half of the selfadjoint matrix C. 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* It corresponds to the level 3 SYRK and level 2 SYR Blas routines. 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath**********************************************************************/ 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs> 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha) 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::conj_if<ConjRhs> cj; 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjLhsType; 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<size; ++i) 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1))) 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez += (alpha * cj(vecY[i])) * ConjLhsType(OtherMap(vecX+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1))); 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs> 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha) 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vecY,vecX,alpha); 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime> 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct selfadjoint_product_selector; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename OtherType, int UpLo> 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true> 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha) 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::blas_traits<OtherType> OtherBlasTraits; 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType; 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType; 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived()); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived()); 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor, 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other; 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(), 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data())); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!UseOtherDirectly) 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo, 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex, 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex> 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualOtherPtr, actualAlpha); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename OtherType, int UpLo> 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false> 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha) 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::blas_traits<OtherType> OtherBlasTraits; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType; 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived()); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived()); 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 }; 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::general_matrix_matrix_triangular_product<Index, 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex, 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex, 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo> 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::run(mat.cols(), actualOther.cols(), 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(), 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.data(), mat.outerStride(), actualAlpha); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// high level API 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename DerivedU> 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez::rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha) 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha); 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SELFADJOINT_PRODUCT_H 124