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