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