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 // 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_SPARSE_SELFADJOINTVIEW_H
00011 #define EIGEN_SPARSE_SELFADJOINTVIEW_H
00012 
00013 namespace Eigen { 
00014 
00029 template<typename Lhs, typename Rhs, int UpLo>
00030 class SparseSelfAdjointTimeDenseProduct;
00031 
00032 template<typename Lhs, typename Rhs, int UpLo>
00033 class DenseTimeSparseSelfAdjointProduct;
00034 
00035 namespace internal {
00036   
00037 template<typename MatrixType, unsigned int UpLo>
00038 struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> {
00039 };
00040 
00041 template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder>
00042 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
00043 
00044 template<int UpLo,typename MatrixType,int DestOrder>
00045 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
00046 
00047 }
00048 
00049 template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
00050   : public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> >
00051 {
00052   public:
00053 
00054     typedef typename MatrixType::Scalar Scalar;
00055     typedef typename MatrixType::Index Index;
00056     typedef Matrix<Index,Dynamic,1> VectorI;
00057     typedef typename MatrixType::Nested MatrixTypeNested;
00058     typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
00059 
00060     inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix)
00061     {
00062       eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices");
00063     }
00064 
00065     inline Index rows() const { return m_matrix.rows(); }
00066     inline Index cols() const { return m_matrix.cols(); }
00067 
00069     const _MatrixTypeNested& matrix() const { return m_matrix; }
00070     _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); }
00071 
00077     template<typename OtherDerived>
00078     SparseSparseProduct<typename OtherDerived::PlainObject, OtherDerived>
00079     operator*(const SparseMatrixBase<OtherDerived>& rhs) const
00080     {
00081       return SparseSparseProduct<typename OtherDerived::PlainObject, OtherDerived>(*this, rhs.derived());
00082     }
00083 
00089     template<typename OtherDerived> friend
00090     SparseSparseProduct<OtherDerived, typename OtherDerived::PlainObject >
00091     operator*(const SparseMatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs)
00092     {
00093       return SparseSparseProduct<OtherDerived, typename OtherDerived::PlainObject>(lhs.derived(), rhs);
00094     }
00095     
00097     template<typename OtherDerived>
00098     SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>
00099     operator*(const MatrixBase<OtherDerived>& rhs) const
00100     {
00101       return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived());
00102     }
00103 
00105     template<typename OtherDerived> friend
00106     DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>
00107     operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs)
00108     {
00109       return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix);
00110     }
00111 
00120     template<typename DerivedU>
00121     SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1));
00122     
00124     template<typename DestScalar,int StorageOrder> void evalTo(SparseMatrix<DestScalar,StorageOrder,Index>& _dest) const
00125     {
00126       internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest);
00127     }
00128     
00129     template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar,ColMajor,Index>& _dest) const
00130     {
00131       // TODO directly evaluate into _dest;
00132       SparseMatrix<DestScalar,ColMajor,Index> tmp(_dest.rows(),_dest.cols());
00133       internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp);
00134       _dest = tmp;
00135     }
00136     
00138     SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const
00139     {
00140       return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm);
00141     }
00142     
00143     template<typename SrcMatrixType,int SrcUpLo>
00144     SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcUpLo>& permutedMatrix)
00145     {
00146       permutedMatrix.evalTo(*this);
00147       return *this;
00148     }
00149 
00150 
00151     SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src)
00152     {
00153       PermutationMatrix<Dynamic> pnull;
00154       return *this = src.twistedBy(pnull);
00155     }
00156 
00157     template<typename SrcMatrixType,unsigned int SrcUpLo>
00158     SparseSelfAdjointView& operator=(const SparseSelfAdjointView<SrcMatrixType,SrcUpLo>& src)
00159     {
00160       PermutationMatrix<Dynamic> pnull;
00161       return *this = src.twistedBy(pnull);
00162     }
00163     
00164 
00165     // const SparseLLT<PlainObject, UpLo> llt() const;
00166     // const SparseLDLT<PlainObject, UpLo> ldlt() const;
00167 
00168   protected:
00169 
00170     typename MatrixType::Nested m_matrix;
00171     mutable VectorI m_countPerRow;
00172     mutable VectorI m_countPerCol;
00173 };
00174 
00175 /***************************************************************************
00176 * Implementation of SparseMatrixBase methods
00177 ***************************************************************************/
00178 
00179 template<typename Derived>
00180 template<unsigned int UpLo>
00181 const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const
00182 {
00183   return derived();
00184 }
00185 
00186 template<typename Derived>
00187 template<unsigned int UpLo>
00188 SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView()
00189 {
00190   return derived();
00191 }
00192 
00193 /***************************************************************************
00194 * Implementation of SparseSelfAdjointView methods
00195 ***************************************************************************/
00196 
00197 template<typename MatrixType, unsigned int UpLo>
00198 template<typename DerivedU>
00199 SparseSelfAdjointView<MatrixType,UpLo>&
00200 SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, const Scalar& alpha)
00201 {
00202   SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint();
00203   if(alpha==Scalar(0))
00204     m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
00205   else
00206     m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>();
00207 
00208   return *this;
00209 }
00210 
00211 /***************************************************************************
00212 * Implementation of sparse self-adjoint time dense matrix
00213 ***************************************************************************/
00214 
00215 namespace internal {
00216 template<typename Lhs, typename Rhs, int UpLo>
00217 struct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> >
00218  : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
00219 {
00220   typedef Dense StorageKind;
00221 };
00222 }
00223 
00224 template<typename Lhs, typename Rhs, int UpLo>
00225 class SparseSelfAdjointTimeDenseProduct
00226   : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
00227 {
00228   public:
00229     EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct)
00230 
00231     SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00232     {}
00233 
00234     template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const
00235     {
00236       EIGEN_ONLY_USED_FOR_DEBUG(alpha);
00237       // TODO use alpha
00238       eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry");
00239       typedef typename internal::remove_all<Lhs>::type _Lhs;
00240       typedef typename _Lhs::InnerIterator LhsInnerIterator;
00241       enum {
00242         LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
00243         ProcessFirstHalf =
00244                  ((UpLo&(Upper|Lower))==(Upper|Lower))
00245               || ( (UpLo&Upper) && !LhsIsRowMajor)
00246               || ( (UpLo&Lower) && LhsIsRowMajor),
00247         ProcessSecondHalf = !ProcessFirstHalf
00248       };
00249       for (Index j=0; j<m_lhs.outerSize(); ++j)
00250       {
00251         LhsInnerIterator i(m_lhs,j);
00252         if (ProcessSecondHalf)
00253         {
00254           while (i && i.index()<j) ++i;
00255           if(i && i.index()==j)
00256           {
00257             dest.row(j) += i.value() * m_rhs.row(j);
00258             ++i;
00259           }
00260         }
00261         for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
00262         {
00263           Index a = LhsIsRowMajor ? j : i.index();
00264           Index b = LhsIsRowMajor ? i.index() : j;
00265           typename Lhs::Scalar v = i.value();
00266           dest.row(a) += (v) * m_rhs.row(b);
00267           dest.row(b) += numext::conj(v) * m_rhs.row(a);
00268         }
00269         if (ProcessFirstHalf && i && (i.index()==j))
00270           dest.row(j) += i.value() * m_rhs.row(j);
00271       }
00272     }
00273 
00274   private:
00275     SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&);
00276 };
00277 
00278 namespace internal {
00279 template<typename Lhs, typename Rhs, int UpLo>
00280 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
00281  : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
00282 {};
00283 }
00284 
00285 template<typename Lhs, typename Rhs, int UpLo>
00286 class DenseTimeSparseSelfAdjointProduct
00287   : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
00288 {
00289   public:
00290     EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct)
00291 
00292     DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00293     {}
00294 
00295     template<typename Dest> void scaleAndAddTo(Dest& /*dest*/, const Scalar& /*alpha*/) const
00296     {
00297       // TODO
00298     }
00299 
00300   private:
00301     DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&);
00302 };
00303 
00304 /***************************************************************************
00305 * Implementation of symmetric copies and permutations
00306 ***************************************************************************/
00307 namespace internal {
00308   
00309 template<typename MatrixType, int UpLo>
00310 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
00311 };
00312 
00313 template<int UpLo,typename MatrixType,int DestOrder>
00314 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
00315 {
00316   typedef typename MatrixType::Index Index;
00317   typedef typename MatrixType::Scalar Scalar;
00318   typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
00319   typedef Matrix<Index,Dynamic,1> VectorI;
00320   
00321   Dest& dest(_dest.derived());
00322   enum {
00323     StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
00324   };
00325   
00326   Index size = mat.rows();
00327   VectorI count;
00328   count.resize(size);
00329   count.setZero();
00330   dest.resize(size,size);
00331   for(Index j = 0; j<size; ++j)
00332   {
00333     Index jp = perm ? perm[j] : j;
00334     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00335     {
00336       Index i = it.index();
00337       Index r = it.row();
00338       Index c = it.col();
00339       Index ip = perm ? perm[i] : i;
00340       if(UpLo==(Upper|Lower))
00341         count[StorageOrderMatch ? jp : ip]++;
00342       else if(r==c)
00343         count[ip]++;
00344       else if(( UpLo==Lower && r>c) || ( UpLo==Upper && r<c))
00345       {
00346         count[ip]++;
00347         count[jp]++;
00348       }
00349     }
00350   }
00351   Index nnz = count.sum();
00352   
00353   // reserve space
00354   dest.resizeNonZeros(nnz);
00355   dest.outerIndexPtr()[0] = 0;
00356   for(Index j=0; j<size; ++j)
00357     dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
00358   for(Index j=0; j<size; ++j)
00359     count[j] = dest.outerIndexPtr()[j];
00360   
00361   // copy data
00362   for(Index j = 0; j<size; ++j)
00363   {
00364     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00365     {
00366       Index i = it.index();
00367       Index r = it.row();
00368       Index c = it.col();
00369       
00370       Index jp = perm ? perm[j] : j;
00371       Index ip = perm ? perm[i] : i;
00372       
00373       if(UpLo==(Upper|Lower))
00374       {
00375         Index k = count[StorageOrderMatch ? jp : ip]++;
00376         dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp;
00377         dest.valuePtr()[k] = it.value();
00378       }
00379       else if(r==c)
00380       {
00381         Index k = count[ip]++;
00382         dest.innerIndexPtr()[k] = ip;
00383         dest.valuePtr()[k] = it.value();
00384       }
00385       else if(( (UpLo&Lower)==Lower && r>c) || ( (UpLo&Upper)==Upper && r<c))
00386       {
00387         if(!StorageOrderMatch)
00388           std::swap(ip,jp);
00389         Index k = count[jp]++;
00390         dest.innerIndexPtr()[k] = ip;
00391         dest.valuePtr()[k] = it.value();
00392         k = count[ip]++;
00393         dest.innerIndexPtr()[k] = jp;
00394         dest.valuePtr()[k] = numext::conj(it.value());
00395       }
00396     }
00397   }
00398 }
00399 
00400 template<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder>
00401 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
00402 {
00403   typedef typename MatrixType::Index Index;
00404   typedef typename MatrixType::Scalar Scalar;
00405   SparseMatrix<Scalar,DstOrder,Index>& dest(_dest.derived());
00406   typedef Matrix<Index,Dynamic,1> VectorI;
00407   enum {
00408     SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor,
00409     StorageOrderMatch = int(SrcOrder) == int(DstOrder),
00410     DstUpLo = DstOrder==RowMajor ? (_DstUpLo==Upper ? Lower : Upper) : _DstUpLo,
00411     SrcUpLo = SrcOrder==RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo
00412   };
00413   
00414   Index size = mat.rows();
00415   VectorI count(size);
00416   count.setZero();
00417   dest.resize(size,size);
00418   for(Index j = 0; j<size; ++j)
00419   {
00420     Index jp = perm ? perm[j] : j;
00421     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00422     {
00423       Index i = it.index();
00424       if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j))
00425         continue;
00426                   
00427       Index ip = perm ? perm[i] : i;
00428       count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
00429     }
00430   }
00431   dest.outerIndexPtr()[0] = 0;
00432   for(Index j=0; j<size; ++j)
00433     dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
00434   dest.resizeNonZeros(dest.outerIndexPtr()[size]);
00435   for(Index j=0; j<size; ++j)
00436     count[j] = dest.outerIndexPtr()[j];
00437   
00438   for(Index j = 0; j<size; ++j)
00439   {
00440     
00441     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00442     {
00443       Index i = it.index();
00444       if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j))
00445         continue;
00446                   
00447       Index jp = perm ? perm[j] : j;
00448       Index ip = perm? perm[i] : i;
00449       
00450       Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
00451       dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp);
00452       
00453       if(!StorageOrderMatch) std::swap(ip,jp);
00454       if( ((int(DstUpLo)==int(Lower) && ip<jp) || (int(DstUpLo)==int(Upper) && ip>jp)))
00455         dest.valuePtr()[k] = numext::conj(it.value());
00456       else
00457         dest.valuePtr()[k] = it.value();
00458     }
00459   }
00460 }
00461 
00462 }
00463 
00464 template<typename MatrixType,int UpLo>
00465 class SparseSymmetricPermutationProduct
00466   : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
00467 {
00468   public:
00469     typedef typename MatrixType::Scalar Scalar;
00470     typedef typename MatrixType::Index Index;
00471   protected:
00472     typedef PermutationMatrix<Dynamic,Dynamic,Index> Perm;
00473   public:
00474     typedef Matrix<Index,Dynamic,1> VectorI;
00475     typedef typename MatrixType::Nested MatrixTypeNested;
00476     typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
00477     
00478     SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm)
00479       : m_matrix(mat), m_perm(perm)
00480     {}
00481     
00482     inline Index rows() const { return m_matrix.rows(); }
00483     inline Index cols() const { return m_matrix.cols(); }
00484     
00485     template<typename DestScalar, int Options, typename DstIndex>
00486     void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest) const
00487     {
00488 //       internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
00489       SparseMatrix<DestScalar,(Options&RowMajor)==RowMajor ? ColMajor : RowMajor, DstIndex> tmp;
00490       internal::permute_symm_to_fullsymm<UpLo>(m_matrix,tmp,m_perm.indices().data());
00491       _dest = tmp;
00492     }
00493     
00494     template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const
00495     {
00496       internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
00497     }
00498     
00499   protected:
00500     MatrixTypeNested m_matrix;
00501     const Perm& m_perm;
00502 
00503 };
00504 
00505 } // end namespace Eigen
00506 
00507 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 12:00:45