SelfadjointProduct.h
Go to the documentation of this file.
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 template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs>
00022 struct selfadjoint_rank1_update;
00023 
00024 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
00025 struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs>
00026 {
00027   static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha)
00028   {
00029     internal::conj_if<ConjRhs> cj;
00030     typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
00031     typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjRhsType;
00032     for (Index i=0; i<size; ++i)
00033     {
00034       Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1)))
00035           += (alpha * cj(vec[i])) * ConjRhsType(OtherMap(vec+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1)));
00036     }
00037   }
00038 };
00039 
00040 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
00041 struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs>
00042 {
00043   static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha)
00044   {
00045     selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vec,alpha);
00046   }
00047 };
00048 
00049 template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime>
00050 struct selfadjoint_product_selector;
00051 
00052 template<typename MatrixType, typename OtherType, int UpLo>
00053 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
00054 {
00055   static void run(MatrixType& mat, const OtherType& other, typename MatrixType::Scalar alpha)
00056   {
00057     typedef typename MatrixType::Scalar Scalar;
00058     typedef typename MatrixType::Index Index;
00059     typedef internal::blas_traits<OtherType> OtherBlasTraits;
00060     typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
00061     typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
00062     typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
00063 
00064     Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
00065 
00066     enum {
00067       StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor,
00068       UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
00069     };
00070     internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
00071 
00072     ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
00073       (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
00074       
00075     if(!UseOtherDirectly)
00076       Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
00077     
00078     selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
00079                               OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
00080                             (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
00081           ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha);
00082   }
00083 };
00084 
00085 template<typename MatrixType, typename OtherType, int UpLo>
00086 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
00087 {
00088   static void run(MatrixType& mat, const OtherType& other, typename MatrixType::Scalar alpha)
00089   {
00090     typedef typename MatrixType::Scalar Scalar;
00091     typedef typename MatrixType::Index Index;
00092     typedef internal::blas_traits<OtherType> OtherBlasTraits;
00093     typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
00094     typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
00095     typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
00096 
00097     Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
00098 
00099     enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
00100 
00101     internal::general_matrix_matrix_triangular_product<Index,
00102       Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor,   OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
00103       Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
00104       MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo>
00105       ::run(mat.cols(), actualOther.cols(),
00106             &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
00107             mat.data(), mat.outerStride(), actualAlpha);
00108   }
00109 };
00110 
00111 // high level API
00112 
00113 template<typename MatrixType, unsigned int UpLo>
00114 template<typename DerivedU>
00115 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
00116 ::rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha)
00117 {
00118   selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
00119 
00120   return *this;
00121 }
00122 
00123 } // end namespace Eigen
00124 
00125 #endif // EIGEN_SELFADJOINT_PRODUCT_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:54