$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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