Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
SelfadjointRank2Update.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_SELFADJOINTRANK2UPTADE_H 00011 #define EIGEN_SELFADJOINTRANK2UPTADE_H 00012 00013 namespace Eigen { 00014 00015 namespace internal { 00016 00017 /* Optimized selfadjoint matrix += alpha * uv' + conj(alpha)*vu' 00018 * It corresponds to the Level2 syr2 BLAS routine 00019 */ 00020 00021 template<typename Scalar, typename Index, typename UType, typename VType, int UpLo> 00022 struct selfadjoint_rank2_update_selector; 00023 00024 template<typename Scalar, typename Index, typename UType, typename VType> 00025 struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Lower> 00026 { 00027 static void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha) 00028 { 00029 const Index size = u.size(); 00030 for (Index i=0; i<size; ++i) 00031 { 00032 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) += 00033 (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.tail(size-i) 00034 + (alpha * numext::conj(v.coeff(i))) * u.tail(size-i); 00035 } 00036 } 00037 }; 00038 00039 template<typename Scalar, typename Index, typename UType, typename VType> 00040 struct selfadjoint_rank2_update_selector<Scalar,Index,UType,VType,Upper> 00041 { 00042 static void run(Scalar* mat, Index stride, const UType& u, const VType& v, const Scalar& alpha) 00043 { 00044 const Index size = u.size(); 00045 for (Index i=0; i<size; ++i) 00046 Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) += 00047 (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.head(i+1) 00048 + (alpha * numext::conj(v.coeff(i))) * u.head(i+1); 00049 } 00050 }; 00051 00052 template<bool Cond, typename T> struct conj_expr_if 00053 : conditional<!Cond, const T&, 00054 CwiseUnaryOp<scalar_conjugate_op<typename traits<T>::Scalar>,T> > {}; 00055 00056 } // end namespace internal 00057 00058 template<typename MatrixType, unsigned int UpLo> 00059 template<typename DerivedU, typename DerivedV> 00060 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> 00061 ::rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha) 00062 { 00063 typedef internal::blas_traits<DerivedU> UBlasTraits; 00064 typedef typename UBlasTraits::DirectLinearAccessType ActualUType; 00065 typedef typename internal::remove_all<ActualUType>::type _ActualUType; 00066 typename internal::add_const_on_value_type<ActualUType>::type actualU = UBlasTraits::extract(u.derived()); 00067 00068 typedef internal::blas_traits<DerivedV> VBlasTraits; 00069 typedef typename VBlasTraits::DirectLinearAccessType ActualVType; 00070 typedef typename internal::remove_all<ActualVType>::type _ActualVType; 00071 typename internal::add_const_on_value_type<ActualVType>::type actualV = VBlasTraits::extract(v.derived()); 00072 00073 // If MatrixType is row major, then we use the routine for lower triangular in the upper triangular case and 00074 // vice versa, and take the complex conjugate of all coefficients and vector entries. 00075 00076 enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 }; 00077 Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived()) 00078 * numext::conj(VBlasTraits::extractScalarFactor(v.derived())); 00079 if (IsRowMajor) 00080 actualAlpha = numext::conj(actualAlpha); 00081 00082 internal::selfadjoint_rank2_update_selector<Scalar, Index, 00083 typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::type>::type, 00084 typename internal::remove_all<typename internal::conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::type>::type, 00085 (IsRowMajor ? int(UpLo==Upper ? Lower : Upper) : UpLo)> 00086 ::run(_expression().const_cast_derived().data(),_expression().outerStride(),actualU,actualV,actualAlpha); 00087 00088 return *this; 00089 } 00090 00091 } // end namespace Eigen 00092 00093 #endif // EIGEN_SELFADJOINTRANK2UPTADE_H
Generated on Tue Jul 12 2022 17:47:00 by 1.7.2