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 
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 
00085     const _MatrixTypeNested& matrix() const { return m_matrix; }
00086     _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); }
00087 
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 
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 
00115     template<typename DerivedU>
00116     SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
00117     
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     
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   //internal::conj_if<SrcUpLo!=DstUpLo> cj;
00371   
00372   Index size = mat.rows();
00373   VectorI count(size);
00374   count.setZero();
00375   dest.resize(size,size);
00376   for(Index j = 0; j<size; ++j)
00377   {
00378     Index jp = perm ? perm[j] : j;
00379     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00380     {
00381       Index i = it.index();
00382       if((SrcUpLo==Lower && i<j) || (SrcUpLo==Upper && i>j))
00383         continue;
00384                   
00385       Index ip = perm ? perm[i] : i;
00386       count[DstUpLo==Lower ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
00387     }
00388   }
00389   dest._outerIndexPtr()[0] = 0;
00390   for(Index j=0; j<size; ++j)
00391     dest._outerIndexPtr()[j+1] = dest._outerIndexPtr()[j] + count[j];
00392   dest.resizeNonZeros(dest._outerIndexPtr()[size]);
00393   for(Index j=0; j<size; ++j)
00394     count[j] = dest._outerIndexPtr()[j];
00395   
00396   for(Index j = 0; j<size; ++j)
00397   {
00398     Index jp = perm ? perm[j] : j;
00399     for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
00400     {
00401       Index i = it.index();
00402       if((SrcUpLo==Lower && i<j) || (SrcUpLo==Upper && i>j))
00403         continue;
00404                   
00405       Index ip = perm? perm[i] : i;
00406       Index k = count[DstUpLo==Lower ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
00407       dest._innerIndexPtr()[k] = DstUpLo==Lower ? (std::max)(ip,jp) : (std::min)(ip,jp);
00408       
00409       if((DstUpLo==Lower && ip<jp) || (DstUpLo==Upper && ip>jp))
00410         dest._valuePtr()[k] = conj(it.value());
00411       else
00412         dest._valuePtr()[k] = it.value();
00413     }
00414   }
00415 }
00416 
00417 }
00418 
00419 template<typename MatrixType,int UpLo>
00420 class SparseSymmetricPermutationProduct
00421   : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
00422 {
00423     typedef PermutationMatrix<Dynamic> Perm;
00424   public:
00425     typedef typename MatrixType::Scalar Scalar;
00426     typedef typename MatrixType::Index Index;
00427     typedef Matrix<Index,Dynamic,1> VectorI;
00428     typedef typename MatrixType::Nested MatrixTypeNested;
00429     typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
00430     
00431     SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm)
00432       : m_matrix(mat), m_perm(perm)
00433     {}
00434     
00435     inline Index rows() const { return m_matrix.rows(); }
00436     inline Index cols() const { return m_matrix.cols(); }
00437     
00438     template<typename DestScalar> void evalTo(SparseMatrix<DestScalar>& _dest) const
00439     {
00440       internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
00441     }
00442     
00443     template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const
00444     {
00445       internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
00446     }
00447     
00448   protected:
00449     const MatrixTypeNested m_matrix;
00450     const Perm& m_perm;
00451 
00452 };
00453 
00454 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:39