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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_SELFADJOINT_PRODUCT_H
00026 #define EIGEN_SELFADJOINT_PRODUCT_H
00027 
00028 /**********************************************************************
00029 * This file implements a self adjoint product: C += A A^T updating only
00030 * half of the selfadjoint matrix C.
00031 * It corresponds to the level 3 SYRK and level 2 SYR Blas routines.
00032 **********************************************************************/
00033 
00034 namespace Eigen { 
00035 
00036 template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs>
00037 struct selfadjoint_rank1_update;
00038 
00039 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
00040 struct selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs>
00041 {
00042   static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha)
00043   {
00044     internal::conj_if<ConjRhs> cj;
00045     typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
00046     typedef typename internal::conditional<ConjLhs,typename OtherMap::ConjugateReturnType,const OtherMap&>::type ConjRhsType;
00047     for (Index i=0; i<size; ++i)
00048     {
00049       Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), (UpLo==Lower ? size-i : (i+1)))
00050           += (alpha * cj(vec[i])) * ConjRhsType(OtherMap(vec+(UpLo==Lower ? i : 0),UpLo==Lower ? size-i : (i+1)));
00051     }
00052   }
00053 };
00054 
00055 template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
00056 struct selfadjoint_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs>
00057 {
00058   static void run(Index size, Scalar* mat, Index stride, const Scalar* vec, Scalar alpha)
00059   {
00060     selfadjoint_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,stride,vec,alpha);
00061   }
00062 };
00063 
00064 template<typename MatrixType, typename OtherType, int UpLo, bool OtherIsVector = OtherType::IsVectorAtCompileTime>
00065 struct selfadjoint_product_selector;
00066 
00067 template<typename MatrixType, typename OtherType, int UpLo>
00068 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
00069 {
00070   static void run(MatrixType& mat, const OtherType& other, typename MatrixType::Scalar alpha)
00071   {
00072     typedef typename MatrixType::Scalar Scalar;
00073     typedef typename MatrixType::Index Index;
00074     typedef internal::blas_traits<OtherType> OtherBlasTraits;
00075     typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
00076     typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
00077     typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
00078 
00079     Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
00080 
00081     enum {
00082       StorageOrder = (internal::traits<MatrixType>::Flags&RowMajorBit) ? RowMajor : ColMajor,
00083       UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
00084     };
00085     internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
00086 
00087     ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
00088       (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
00089       
00090     if(!UseOtherDirectly)
00091       Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
00092     
00093     selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
00094                               OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
00095                             (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
00096           ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha);
00097   }
00098 };
00099 
00100 template<typename MatrixType, typename OtherType, int UpLo>
00101 struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,false>
00102 {
00103   static void run(MatrixType& mat, const OtherType& other, typename MatrixType::Scalar alpha)
00104   {
00105     typedef typename MatrixType::Scalar Scalar;
00106     typedef typename MatrixType::Index Index;
00107     typedef internal::blas_traits<OtherType> OtherBlasTraits;
00108     typedef typename OtherBlasTraits::DirectLinearAccessType ActualOtherType;
00109     typedef typename internal::remove_all<ActualOtherType>::type _ActualOtherType;
00110     typename internal::add_const_on_value_type<ActualOtherType>::type actualOther = OtherBlasTraits::extract(other.derived());
00111 
00112     Scalar actualAlpha = alpha * OtherBlasTraits::extractScalarFactor(other.derived());
00113 
00114     enum { IsRowMajor = (internal::traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
00115 
00116     internal::general_matrix_matrix_triangular_product<Index,
00117       Scalar, _ActualOtherType::Flags&RowMajorBit ? RowMajor : ColMajor,   OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
00118       Scalar, _ActualOtherType::Flags&RowMajorBit ? ColMajor : RowMajor, (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex,
00119       MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo>
00120       ::run(mat.cols(), actualOther.cols(),
00121             &actualOther.coeffRef(0,0), actualOther.outerStride(), &actualOther.coeffRef(0,0), actualOther.outerStride(),
00122             mat.data(), mat.outerStride(), actualAlpha);
00123   }
00124 };
00125 
00126 // high level API
00127 
00128 template<typename MatrixType, unsigned int UpLo>
00129 template<typename DerivedU>
00130 SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
00131 ::rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha)
00132 {
00133   selfadjoint_product_selector<MatrixType,DerivedU,UpLo>::run(_expression().const_cast_derived(), u.derived(), alpha);
00134 
00135   return *this;
00136 }
00137 
00138 } // end namespace Eigen
00139 
00140 #endif // EIGEN_SELFADJOINT_PRODUCT_H