1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5// 6// This Source Code Form is subject to the terms of the Mozilla 7// Public License v. 2.0. If a copy of the MPL was not distributed 8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10#ifndef EIGEN_SPARSE_SELFADJOINTVIEW_H 11#define EIGEN_SPARSE_SELFADJOINTVIEW_H 12 13namespace Eigen { 14 15/** \ingroup SparseCore_Module 16 * \class SparseSelfAdjointView 17 * 18 * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix. 19 * 20 * \param MatrixType the type of the dense matrix storing the coefficients 21 * \param UpLo can be either \c #Lower or \c #Upper 22 * 23 * This class is an expression of a sefladjoint matrix from a triangular part of a matrix 24 * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() 25 * and most of the time this is the only way that it is used. 26 * 27 * \sa SparseMatrixBase::selfadjointView() 28 */ 29template<typename Lhs, typename Rhs, int UpLo> 30class SparseSelfAdjointTimeDenseProduct; 31 32template<typename Lhs, typename Rhs, int UpLo> 33class DenseTimeSparseSelfAdjointProduct; 34 35namespace internal { 36 37template<typename MatrixType, unsigned int UpLo> 38struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> { 39}; 40 41template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder> 42void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0); 43 44template<int UpLo,typename MatrixType,int DestOrder> 45void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0); 46 47} 48 49template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView 50 : public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> > 51{ 52 public: 53 54 typedef typename MatrixType::Scalar Scalar; 55 typedef typename MatrixType::Index Index; 56 typedef Matrix<Index,Dynamic,1> VectorI; 57 typedef typename MatrixType::Nested MatrixTypeNested; 58 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested; 59 60 inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix) 61 { 62 eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices"); 63 } 64 65 inline Index rows() const { return m_matrix.rows(); } 66 inline Index cols() const { return m_matrix.cols(); } 67 68 /** \internal \returns a reference to the nested matrix */ 69 const _MatrixTypeNested& matrix() const { return m_matrix; } 70 _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); } 71 72 /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ 73 template<typename OtherDerived> 74 SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo> 75 operator*(const MatrixBase<OtherDerived>& rhs) const 76 { 77 return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived()); 78 } 79 80 /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ 81 template<typename OtherDerived> friend 82 DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo> 83 operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs) 84 { 85 return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix); 86 } 87 88 /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: 89 * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. 90 * 91 * \returns a reference to \c *this 92 * 93 * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply 94 * call this function with u.adjoint(). 95 */ 96 template<typename DerivedU> 97 SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1)); 98 99 /** \internal triggered by sparse_matrix = SparseSelfadjointView; */ 100 template<typename DestScalar,int StorageOrder> void evalTo(SparseMatrix<DestScalar,StorageOrder,Index>& _dest) const 101 { 102 internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest); 103 } 104 105 template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar,ColMajor,Index>& _dest) const 106 { 107 // TODO directly evaluate into _dest; 108 SparseMatrix<DestScalar,ColMajor,Index> tmp(_dest.rows(),_dest.cols()); 109 internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp); 110 _dest = tmp; 111 } 112 113 /** \returns an expression of P H P^-1 */ 114 SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const 115 { 116 return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm); 117 } 118 119 template<typename SrcMatrixType,int SrcUpLo> 120 SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcUpLo>& permutedMatrix) 121 { 122 permutedMatrix.evalTo(*this); 123 return *this; 124 } 125 126 127 SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) 128 { 129 PermutationMatrix<Dynamic> pnull; 130 return *this = src.twistedBy(pnull); 131 } 132 133 template<typename SrcMatrixType,unsigned int SrcUpLo> 134 SparseSelfAdjointView& operator=(const SparseSelfAdjointView<SrcMatrixType,SrcUpLo>& src) 135 { 136 PermutationMatrix<Dynamic> pnull; 137 return *this = src.twistedBy(pnull); 138 } 139 140 141 // const SparseLLT<PlainObject, UpLo> llt() const; 142 // const SparseLDLT<PlainObject, UpLo> ldlt() const; 143 144 protected: 145 146 typename MatrixType::Nested m_matrix; 147 mutable VectorI m_countPerRow; 148 mutable VectorI m_countPerCol; 149}; 150 151/*************************************************************************** 152* Implementation of SparseMatrixBase methods 153***************************************************************************/ 154 155template<typename Derived> 156template<unsigned int UpLo> 157const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const 158{ 159 return derived(); 160} 161 162template<typename Derived> 163template<unsigned int UpLo> 164SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() 165{ 166 return derived(); 167} 168 169/*************************************************************************** 170* Implementation of SparseSelfAdjointView methods 171***************************************************************************/ 172 173template<typename MatrixType, unsigned int UpLo> 174template<typename DerivedU> 175SparseSelfAdjointView<MatrixType,UpLo>& 176SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha) 177{ 178 SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint(); 179 if(alpha==Scalar(0)) 180 m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>(); 181 else 182 m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>(); 183 184 return *this; 185} 186 187/*************************************************************************** 188* Implementation of sparse self-adjoint time dense matrix 189***************************************************************************/ 190 191namespace internal { 192template<typename Lhs, typename Rhs, int UpLo> 193struct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> > 194 : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> > 195{ 196 typedef Dense StorageKind; 197}; 198} 199 200template<typename Lhs, typename Rhs, int UpLo> 201class SparseSelfAdjointTimeDenseProduct 202 : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> 203{ 204 public: 205 EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct) 206 207 SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) 208 {} 209 210 template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const 211 { 212 // TODO use alpha 213 eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry"); 214 typedef typename internal::remove_all<Lhs>::type _Lhs; 215 typedef typename internal::remove_all<Rhs>::type _Rhs; 216 typedef typename _Lhs::InnerIterator LhsInnerIterator; 217 enum { 218 LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, 219 ProcessFirstHalf = 220 ((UpLo&(Upper|Lower))==(Upper|Lower)) 221 || ( (UpLo&Upper) && !LhsIsRowMajor) 222 || ( (UpLo&Lower) && LhsIsRowMajor), 223 ProcessSecondHalf = !ProcessFirstHalf 224 }; 225 for (Index j=0; j<m_lhs.outerSize(); ++j) 226 { 227 LhsInnerIterator i(m_lhs,j); 228 if (ProcessSecondHalf) 229 { 230 while (i && i.index()<j) ++i; 231 if(i && i.index()==j) 232 { 233 dest.row(j) += i.value() * m_rhs.row(j); 234 ++i; 235 } 236 } 237 for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) 238 { 239 Index a = LhsIsRowMajor ? j : i.index(); 240 Index b = LhsIsRowMajor ? i.index() : j; 241 typename Lhs::Scalar v = i.value(); 242 dest.row(a) += (v) * m_rhs.row(b); 243 dest.row(b) += internal::conj(v) * m_rhs.row(a); 244 } 245 if (ProcessFirstHalf && i && (i.index()==j)) 246 dest.row(j) += i.value() * m_rhs.row(j); 247 } 248 } 249 250 private: 251 SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&); 252}; 253 254namespace internal { 255template<typename Lhs, typename Rhs, int UpLo> 256struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> > 257 : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> > 258{}; 259} 260 261template<typename Lhs, typename Rhs, int UpLo> 262class DenseTimeSparseSelfAdjointProduct 263 : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> 264{ 265 public: 266 EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct) 267 268 DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) 269 {} 270 271 template<typename Dest> void scaleAndAddTo(Dest& /*dest*/, Scalar /*alpha*/) const 272 { 273 // TODO 274 } 275 276 private: 277 DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&); 278}; 279 280/*************************************************************************** 281* Implementation of symmetric copies and permutations 282***************************************************************************/ 283namespace internal { 284 285template<typename MatrixType, int UpLo> 286struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> { 287}; 288 289template<int UpLo,typename MatrixType,int DestOrder> 290void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm) 291{ 292 typedef typename MatrixType::Index Index; 293 typedef typename MatrixType::Scalar Scalar; 294 typedef SparseMatrix<Scalar,DestOrder,Index> Dest; 295 typedef Matrix<Index,Dynamic,1> VectorI; 296 297 Dest& dest(_dest.derived()); 298 enum { 299 StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor) 300 }; 301 302 Index size = mat.rows(); 303 VectorI count; 304 count.resize(size); 305 count.setZero(); 306 dest.resize(size,size); 307 for(Index j = 0; j<size; ++j) 308 { 309 Index jp = perm ? perm[j] : j; 310 for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 311 { 312 Index i = it.index(); 313 Index r = it.row(); 314 Index c = it.col(); 315 Index ip = perm ? perm[i] : i; 316 if(UpLo==(Upper|Lower)) 317 count[StorageOrderMatch ? jp : ip]++; 318 else if(r==c) 319 count[ip]++; 320 else if(( UpLo==Lower && r>c) || ( UpLo==Upper && r<c)) 321 { 322 count[ip]++; 323 count[jp]++; 324 } 325 } 326 } 327 Index nnz = count.sum(); 328 329 // reserve space 330 dest.resizeNonZeros(nnz); 331 dest.outerIndexPtr()[0] = 0; 332 for(Index j=0; j<size; ++j) 333 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j]; 334 for(Index j=0; j<size; ++j) 335 count[j] = dest.outerIndexPtr()[j]; 336 337 // copy data 338 for(Index j = 0; j<size; ++j) 339 { 340 for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 341 { 342 Index i = it.index(); 343 Index r = it.row(); 344 Index c = it.col(); 345 346 Index jp = perm ? perm[j] : j; 347 Index ip = perm ? perm[i] : i; 348 349 if(UpLo==(Upper|Lower)) 350 { 351 Index k = count[StorageOrderMatch ? jp : ip]++; 352 dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp; 353 dest.valuePtr()[k] = it.value(); 354 } 355 else if(r==c) 356 { 357 Index k = count[ip]++; 358 dest.innerIndexPtr()[k] = ip; 359 dest.valuePtr()[k] = it.value(); 360 } 361 else if(( (UpLo&Lower)==Lower && r>c) || ( (UpLo&Upper)==Upper && r<c)) 362 { 363 if(!StorageOrderMatch) 364 std::swap(ip,jp); 365 Index k = count[jp]++; 366 dest.innerIndexPtr()[k] = ip; 367 dest.valuePtr()[k] = it.value(); 368 k = count[ip]++; 369 dest.innerIndexPtr()[k] = jp; 370 dest.valuePtr()[k] = internal::conj(it.value()); 371 } 372 } 373 } 374} 375 376template<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder> 377void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm) 378{ 379 typedef typename MatrixType::Index Index; 380 typedef typename MatrixType::Scalar Scalar; 381 SparseMatrix<Scalar,DstOrder,Index>& dest(_dest.derived()); 382 typedef Matrix<Index,Dynamic,1> VectorI; 383 enum { 384 SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor, 385 StorageOrderMatch = int(SrcOrder) == int(DstOrder), 386 DstUpLo = DstOrder==RowMajor ? (_DstUpLo==Upper ? Lower : Upper) : _DstUpLo, 387 SrcUpLo = SrcOrder==RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo 388 }; 389 390 Index size = mat.rows(); 391 VectorI count(size); 392 count.setZero(); 393 dest.resize(size,size); 394 for(Index j = 0; j<size; ++j) 395 { 396 Index jp = perm ? perm[j] : j; 397 for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 398 { 399 Index i = it.index(); 400 if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j)) 401 continue; 402 403 Index ip = perm ? perm[i] : i; 404 count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; 405 } 406 } 407 dest.outerIndexPtr()[0] = 0; 408 for(Index j=0; j<size; ++j) 409 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j]; 410 dest.resizeNonZeros(dest.outerIndexPtr()[size]); 411 for(Index j=0; j<size; ++j) 412 count[j] = dest.outerIndexPtr()[j]; 413 414 for(Index j = 0; j<size; ++j) 415 { 416 417 for(typename MatrixType::InnerIterator it(mat,j); it; ++it) 418 { 419 Index i = it.index(); 420 if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j)) 421 continue; 422 423 Index jp = perm ? perm[j] : j; 424 Index ip = perm? perm[i] : i; 425 426 Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; 427 dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp); 428 429 if(!StorageOrderMatch) std::swap(ip,jp); 430 if( ((int(DstUpLo)==int(Lower) && ip<jp) || (int(DstUpLo)==int(Upper) && ip>jp))) 431 dest.valuePtr()[k] = conj(it.value()); 432 else 433 dest.valuePtr()[k] = it.value(); 434 } 435 } 436} 437 438} 439 440template<typename MatrixType,int UpLo> 441class SparseSymmetricPermutationProduct 442 : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> > 443{ 444 public: 445 typedef typename MatrixType::Scalar Scalar; 446 typedef typename MatrixType::Index Index; 447 protected: 448 typedef PermutationMatrix<Dynamic,Dynamic,Index> Perm; 449 public: 450 typedef Matrix<Index,Dynamic,1> VectorI; 451 typedef typename MatrixType::Nested MatrixTypeNested; 452 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested; 453 454 SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm) 455 : m_matrix(mat), m_perm(perm) 456 {} 457 458 inline Index rows() const { return m_matrix.rows(); } 459 inline Index cols() const { return m_matrix.cols(); } 460 461 template<typename DestScalar, int Options, typename DstIndex> 462 void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest) const 463 { 464 internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data()); 465 } 466 467 template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const 468 { 469 internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data()); 470 } 471 472 protected: 473 MatrixTypeNested m_matrix; 474 const Perm& m_perm; 475 476}; 477 478} // end namespace Eigen 479 480#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H 481