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