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


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:26:20