Main MRPT website > C++ reference
MRPT logo

SparseSelfAdjointView.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_SPARSE_SELFADJOINTVIEW_H
00026 #define EIGEN_SPARSE_SELFADJOINTVIEW_H
00027 
00028 /** \class SparseSelfAdjointView
00029   *
00030   *
00031   * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
00032   *
00033   * \param MatrixType the type of the dense matrix storing the coefficients
00034   * \param UpLo can be either \c Lower or \c Upper
00035   *
00036   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
00037   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
00038   * and most of the time this is the only way that it is used.
00039   *
00040   * \sa SparseMatrixBase::selfAdjointView()
00041   */
00042 template<typename Lhs, typename Rhs, int UpLo>
00043 class SparseSelfAdjointTimeDenseProduct;
00044 
00045 template<typename Lhs, typename Rhs, int UpLo>
00046 class DenseTimeSparseSelfAdjointProduct;
00047 
00048 template<typename MatrixType,int UpLo>
00049 class SparseSymmetricPermutationProduct;
00050 
00051 namespace internal {
00052   
00053 template<typename MatrixType, unsigned int UpLo>
00054 struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> {
00055 };
00056 
00057 template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder>
00058 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
00059 
00060 template<int UpLo,typename MatrixType,int DestOrder>
00061 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
00062 
00063 }
00064 
00065 template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
00066   : public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> >
00067 {
00068   public:
00069 
00070     typedef typename MatrixType::Scalar Scalar;
00071     typedef typename MatrixType::Index Index;
00072     typedef Matrix<Index,Dynamic,1> VectorI;
00073     typedef typename MatrixType::Nested MatrixTypeNested;
00074     typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
00075 
00076     inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix)
00077     {
00078       eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices");
00079     }
00080 
00081     inline Index rows() const { return m_matrix.rows(); }
00082     inline Index cols() const { return m_matrix.cols(); }
00083 
00084     /** \internal \returns a reference to the nested matrix */
00085     const _MatrixTypeNested& matrix() const { return m_matrix; }
00086     _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); }
00087 
00088     /** Efficient sparse self-adjoint matrix times dense vector/matrix product */
00089     template<typename OtherDerived>
00090     SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>
00091     operator*(const MatrixBase<OtherDerived>& rhs) const
00092     {
00093       return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived());
00094     }
00095 
00096     /** Efficient dense vector/matrix times sparse self-adjoint matrix product */
00097     template<typename OtherDerived> friend
00098     DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>
00099     operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs)
00100     {
00101       return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix);
00102     }
00103 
00104     /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
00105       * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
00106       *
00107       * \returns a reference to \c *this
00108       *
00109       * Note that it is faster to set alpha=0 than initializing the matrix to zero
00110       * and then keep the default value alpha=1.
00111       *
00112       * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
00113       * call this function with u.adjoint().
00114       */
00115     template<typename DerivedU>
00116     SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
00117     
00118     /** \internal triggered by sparse_matrix = SparseSelfadjointView; */
00119     template<typename DestScalar> void evalTo(SparseMatrix<DestScalar>& _dest) const
00120     {
00121       internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest);
00122     }
00123     
00124     template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar>& _dest) const
00125     {
00126       // TODO directly evaluate into _dest;
00127       SparseMatrix<DestScalar> tmp(_dest.rows(),_dest.cols());
00128       internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp);
00129       _dest = tmp;
00130     }
00131     
00132     /** \returns an expression of P^-1 H P */
00133     SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix<Dynamic>& perm) const
00134     {
00135       return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm);
00136     }
00137     
00138     template<typename SrcMatrixType,int SrcUpLo>
00139     SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcUpLo>& permutedMatrix)
00140     {
00141       permutedMatrix.evalTo(*this);
00142       return *this;
00143     }
00144     
00145 
00146     // const SparseLLT<PlainObject, UpLo> llt() const;
00147     // const SparseLDLT<PlainObject, UpLo> ldlt() const;
00148 
00149   protected:
00150 
00151     const typename MatrixType::Nested m_matrix;
00152     mutable VectorI m_countPerRow;
00153     mutable VectorI m_countPerCol;
00154 };
00155 
00156 /***************************************************************************
00157 * Implementation of SparseMatrixBase methods
00158 ***************************************************************************/
00159 
00160 template<typename Derived>
00161 template<unsigned int UpLo>
00162 const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const
00163 {
00164   return derived();
00165 }
00166 
00167 template<typename Derived>
00168 template<unsigned int UpLo>
00169 SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView()
00170 {
00171   return derived();
00172 }
00173 
00174 /***************************************************************************
00175 * Implementation of SparseSelfAdjointView methods
00176 ***************************************************************************/
00177 
00178 template<typename MatrixType, unsigned int UpLo>
00179 template<typename DerivedU>
00180 SparseSelfAdjointView<MatrixType,UpLo>&
00181 SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha)
00182 {
00183   SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint();
00184   if(alpha==Scalar(0))
00185     m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
00186   else
00187     m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>();
00188 
00189   return *this;
00190 }
00191 
00192 /***************************************************************************
00193 * Implementation of sparse self-adjoint time dense matrix
00194 ***************************************************************************/
00195 
00196 namespace internal {
00197 template<typename Lhs, typename Rhs, int UpLo>
00198 struct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> >
00199  : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
00200 {
00201   typedef Dense StorageKind;
00202 };
00203 }
00204 
00205 template<typename Lhs, typename Rhs, int UpLo>
00206 class SparseSelfAdjointTimeDenseProduct
00207   : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
00208 {
00209   public:
00210     EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct)
00211 
00212     SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00213     {}
00214 
00215     template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
00216     {
00217       // TODO use alpha
00218       eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry");
00219       typedef typename internal::remove_all<Lhs>::type _Lhs;
00220       typedef typename internal::remove_all<Rhs>::type _Rhs;
00221       typedef typename _Lhs::InnerIterator LhsInnerIterator;
00222       enum {
00223         LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
00224         ProcessFirstHalf =
00225                  ((UpLo&(Upper|Lower))==(Upper|Lower))
00226               || ( (UpLo&Upper) && !LhsIsRowMajor)
00227               || ( (UpLo&Lower) && LhsIsRowMajor),
00228         ProcessSecondHalf = !ProcessFirstHalf
00229       };
00230       for (Index j=0; j<m_lhs.outerSize(); ++j)
00231       {
00232         LhsInnerIterator i(m_lhs,j);
00233         if (ProcessSecondHalf && i && (i.index()==j))
00234         {
00235           dest.row(j) += i.value() * m_rhs.row(j);
00236           ++i;
00237         }
00238         Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0));
00239         for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
00240         {
00241           Index a = LhsIsRowMajor ? j : i.index();
00242           Index b = LhsIsRowMajor ? i.index() : j;
00243           typename Lhs::Scalar v = i.value();
00244           dest.row(a) += (v) * m_rhs.row(b);
00245           dest.row(b) += internal::conj(v) * m_rhs.row(a);
00246         }
00247         if (ProcessFirstHalf && i && (i.index()==j))
00248           dest.row(j) += i.value() * m_rhs.row(j);
00249       }
00250     }
00251 
00252   private:
00253     SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&);
00254 };
00255 
00256 namespace internal {
00257 template<typename Lhs, typename Rhs, int UpLo>
00258 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
00259  : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
00260 {};
00261 }
00262 
00263 template<typename Lhs, typename Rhs, int UpLo>
00264 class DenseTimeSparseSelfAdjointProduct
00265   : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
00266 {
00267   public:
00268     EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct)
00269 
00270     DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00271     {}
00272 
00273     template<typename Dest> void scaleAndAddTo(Dest& /*dest*/, Scalar /*alpha*/) const
00274     {
00275       // TODO
00276     }
00277 
00278   private:
00279     DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&);
00280 };
00281 
00282 /***************************************************************************
00283 * Implementation of symmetric copies and permutations
00284 ***************************************************************************/
00285 namespace internal {
00286   
00287 template<typename MatrixType, int UpLo>
00288 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
00289 };
00290 
00291 template<int UpLo,typename MatrixType,int DestOrder>
00292 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
00293 {
00294   typedef typename MatrixType::Index Index;
00295   typedef typename MatrixType::Scalar Scalar;
00296   typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
00297   typedef Matrix<Index,Dynamic,1> VectorI;
00298   
00299   Dest& dest(_dest.derived());
00300   enum {
00301     StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
00302   };
00303   eigen_assert(perm==0);
00304   Index size = mat.rows();
00305   VectorI count;
00306   count.resize(size);
00307   count.setZero();
00308   dest.resize(size,size);
00309   for(Index j = 0; j<size; ++j)
00310   {
00311     Index jp = perm ? perm[j] : j;
00312     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00313     {
00314       Index i = it.index();
00315       Index ip = perm ? perm[i] : i;
00316       if(i==j)
00317         count[ip]++;
00318       else if((UpLo==Lower && i>j) || (UpLo==Upper && i<j))
00319       {
00320         count[ip]++;
00321         count[jp]++;
00322       }
00323     }
00324   }
00325   Index nnz = count.sum();
00326   
00327   // reserve space
00328   dest.reserve(nnz);
00329   dest._outerIndexPtr()[0] = 0;
00330   for(Index j=0; j<size; ++j)
00331     dest._outerIndexPtr()[j+1] = dest._outerIndexPtr()[j] + count[j];
00332   for(Index j=0; j<size; ++j)
00333     count[j] = dest._outerIndexPtr()[j];
00334   
00335   // copy data
00336   for(Index j = 0; j<size; ++j)
00337   {
00338     Index jp = perm ? perm[j] : j;
00339     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00340     {
00341       Index i = it.index();
00342       Index ip = perm ? perm[i] : i;
00343       if(i==j)
00344       {
00345         int k = count[ip]++;
00346         dest._innerIndexPtr()[k] = ip;
00347         dest._valuePtr()[k] = it.value();
00348       }
00349       else if((UpLo==Lower && i>j) || (UpLo==Upper && i<j))
00350       {
00351         int k = count[jp]++;
00352         dest._innerIndexPtr()[k] = ip;
00353         dest._valuePtr()[k] = it.value();
00354         k = count[ip]++;
00355         dest._innerIndexPtr()[k] = jp;
00356         dest._valuePtr()[k] = internal::conj(it.value());
00357       }
00358     }
00359   }
00360 }
00361 
00362 template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder>
00363 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
00364 {
00365   typedef typename MatrixType::Index Index;
00366   typedef typename MatrixType::Scalar Scalar;
00367   typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
00368   Dest& dest(_dest.derived());
00369   typedef Matrix<Index,Dynamic,1> VectorI;
00370   
00371   Index size = mat.rows();
00372   VectorI count(size);
00373   count.setZero();
00374   dest.resize(size,size);
00375   for(Index j = 0; j<size; ++j)
00376   {
00377     Index jp = perm ? perm[j] : j;
00378     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00379     {
00380       Index i = it.index();
00381       if((SrcUpLo==Lower && i<j) || (SrcUpLo==Upper && i>j))
00382         continue;
00383                   
00384       Index ip = perm ? perm[i] : i;
00385       count[DstUpLo==Lower ? std::min(ip,jp) : std::max(ip,jp)]++;
00386     }
00387   }
00388   dest._outerIndexPtr()[0] = 0;
00389   for(Index j=0; j<size; ++j)
00390     dest._outerIndexPtr()[j+1] = dest._outerIndexPtr()[j] + count[j];
00391   dest.resizeNonZeros(dest._outerIndexPtr()[size]);
00392   for(Index j=0; j<size; ++j)
00393     count[j] = dest._outerIndexPtr()[j];
00394   
00395   for(Index j = 0; j<size; ++j)
00396   {
00397     Index jp = perm ? perm[j] : j;
00398     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00399     {
00400       Index i = it.index();
00401       if((SrcUpLo==Lower && i<j) || (SrcUpLo==Upper && i>j))
00402         continue;
00403                   
00404       Index ip = perm? perm[i] : i;
00405       Index k = count[DstUpLo==Lower ? std::min(ip,jp) : std::max(ip,jp)]++;
00406       dest._innerIndexPtr()[k] = DstUpLo==Lower ? std::max(ip,jp) : std::min(ip,jp);
00407       dest._valuePtr()[k] = it.value();
00408     }
00409   }
00410 }
00411 
00412 }
00413 
00414 template<typename MatrixType,int UpLo>
00415 class SparseSymmetricPermutationProduct
00416   : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
00417 {
00418     typedef PermutationMatrix<Dynamic> Perm;
00419   public:
00420     typedef typename MatrixType::Scalar Scalar;
00421     typedef typename MatrixType::Index Index;
00422     typedef Matrix<Index,Dynamic,1> VectorI;
00423     typedef typename MatrixType::Nested MatrixTypeNested;
00424     typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
00425     
00426     SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm)
00427       : m_matrix(mat), m_perm(perm)
00428     {}
00429     
00430     inline Index rows() const { return m_matrix.rows(); }
00431     inline Index cols() const { return m_matrix.cols(); }
00432     
00433     template<typename DestScalar> void evalTo(SparseMatrix<DestScalar>& _dest) const
00434     {
00435       internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
00436     }
00437     
00438     template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const
00439     {
00440       internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
00441     }
00442     
00443   protected:
00444     const MatrixTypeNested m_matrix;
00445     const Perm& m_perm;
00446 
00447 };
00448 
00449 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:40:17 UTC 2011