Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
SelfadjointProduct.h
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_SELFADJOINT_PRODUCT_H 00011 #define EIGEN_SELFADJOINT_PRODUCT_H 00012 00013 /********************************************************************** 00014 * This file implements a self adjoint product: C += A A^T updating only 00015 * half of the selfadjoint matrix C. 00016 * It corresponds to the level 3 SYRK and level 2 SYR Blas routines. 00017 **********************************************************************/ 00018 00019 namespace Eigen { 00020 00021 00022 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> 00023 struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs> 00024 { 00025 static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha) 00026 { 00027 internal::conj_if<ConjRhs> cj; 00028 typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap; 00029 typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjLhsType; 00030 for (Index i=0; i<size; ++i) 00031 { 00032 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1))) 00033 += (alpha * cj(vecY[i])) * ConjLhsType(OtherMap(vecX+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1))); 00034 } 00035 } 00036 }; 00037 00038 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs> 00039 struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs> 00040 { 00041 static void run(Index size, Scalar* mat, Index stride, const Scalar* vecX, const Scalar* vecY, const Scalar& alpha) 00042 { 00043 selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vecY,vecX,alpha); 00044 } 00045 }; 00046 00047 template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime> 00048 struct selfadjoint_product_selector; 00049 00050 template<typename MatrixType, typename OtherType, int UpLo> 00051 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true> 00052 { 00053 static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha) 00054 { 00055 typedef typename MatrixType::Scalar Scalar; 00056 typedef typename MatrixType::Index Index; 00057 typedef internal::blas_traits<OtherType> OtherBlasTraits; 00058 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType; 00059 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType; 00060 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived()); 00061 00062 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived()); 00063 00064 enum { 00065 StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor, 00066 UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1 00067 }; 00068 internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other; 00069 00070 ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(), 00071 (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data())); 00072 00073 if(!UseOtherDirectly) 00074 Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther; 00075 00076 selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo, 00077 OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex, 00078 (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex> 00079 ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualOtherPtr, actualAlpha); 00080 } 00081 }; 00082 00083 template<typename MatrixType, typename OtherType, int UpLo> 00084 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false> 00085 { 00086 static void run(MatrixType& mat, const OtherType& other, const typename MatrixType::Scalar& alpha) 00087 { 00088 typedef typename MatrixType::Scalar Scalar; 00089 typedef typename MatrixType::Index Index; 00090 typedef internal::blas_traits<OtherType> OtherBlasTraits; 00091 typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType; 00092 typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType; 00093 typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived()); 00094 00095 Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived()); 00096 00097 enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 }; 00098 00099 internal::general_matrix_matrix_triangular_product<Index, 00100 Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor, OtherBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex, 00101 Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex, 00102 MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo> 00103 ::run(mat.cols(), actualOther.cols(), 00104 &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(), 00105 mat.data(), mat.outerStride(), actualAlpha); 00106 } 00107 }; 00108 00109 // high level API 00110 00111 template<typename MatrixType, unsigned int UpLo> 00112 template<typename DerivedU> 00113 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> 00114 ::rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha) 00115 { 00116 selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha); 00117 00118 return *this; 00119 } 00120 00121 } // end namespace Eigen 00122 00123 #endif // EIGEN_SELFADJOINT_PRODUCT_H
Generated on Tue Jul 12 2022 17:47:00 by
![doxygen](doxygen.png)