SparseSelfAdjointView.h
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
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
69 const _MatrixTypeNested& matrix() const { return m_matrix; }
70 _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); }
71
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
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
96 template<typename DerivedU>
97 SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
98
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
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 EIGEN_ONLY_USED_FOR_DEBUG(alpha);
213 // TODO use alpha
214 eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry");
215 typedef typename internal::remove_all<Lhs>::type _Lhs;
216 typedef typename internal::remove_all<Rhs>::type _Rhs;
217 typedef typename _Lhs::InnerIterator LhsInnerIterator;
218 enum {
219 LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
220 ProcessFirstHalf =
221 ((UpLo&(Upper|Lower))==(Upper|Lower))
222 || ( (UpLo&Upper) && !LhsIsRowMajor)
223 || ( (UpLo&Lower) && LhsIsRowMajor),
224 ProcessSecondHalf = !ProcessFirstHalf
225 };
226 for (Index j=0; j<m_lhs.outerSize(); ++j)
227 {
228 LhsInnerIterator i(m_lhs,j);
229 if (ProcessSecondHalf)
230 {
231 while (i && i.index()<j) ++i;
232 if(i && i.index()==j)
233 {
234 dest.row(j) += i.value() * m_rhs.row(j);
235 ++i;
236 }
237 }
238 for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
239 {
240 Index a = LhsIsRowMajor ? j : i.index();
241 Index b = LhsIsRowMajor ? i.index() : j;
242 typename Lhs::Scalar v = i.value();
243 dest.row(a) += (v) * m_rhs.row(b);
244 dest.row(b) += internal::conj(v) * m_rhs.row(a);
245 }
246 if (ProcessFirstHalf && i && (i.index()==j))
247 dest.row(j) += i.value() * m_rhs.row(j);
248 }
249 }
250
251 private:
252 SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&);
253};
254
255namespace internal {
256template<typename Lhs, typename Rhs, int UpLo>
257struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
258 : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
259{};
260}
261
262template<typename Lhs, typename Rhs, int UpLo>
263class DenseTimeSparseSelfAdjointProduct
264 : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
265{
266 public:
267 EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct)
268
269 DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
270 {}
271
272 template<typename Dest> void scaleAndAddTo(Dest& /*dest*/, Scalar /*alpha*/) const
273 {
274 // TODO
275 }
276
277 private:
278 DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&);
279};
280
281/***************************************************************************
282* Implementation of symmetric copies and permutations
283***************************************************************************/
284namespace internal {
285
286template<typename MatrixType, int UpLo>
287struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
288};
289
290template<int UpLo,typename MatrixType,int DestOrder>
291void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
292{
293 typedef typename MatrixType::Index Index;
294 typedef typename MatrixType::Scalar Scalar;
295 typedef SparseMatrix<Scalar,DestOrder,Index> Dest;
296 typedef Matrix<Index,Dynamic,1> VectorI;
297
298 Dest& dest(_dest.derived());
299 enum {
300 StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
301 };
302
303 Index size = mat.rows();
304 VectorI count;
305 count.resize(size);
306 count.setZero();
307 dest.resize(size,size);
308 for(Index j = 0; j<size; ++j)
309 {
310 Index jp = perm ? perm[j] : j;
311 for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
312 {
313 Index i = it.index();
314 Index r = it.row();
315 Index c = it.col();
316 Index ip = perm ? perm[i] : i;
317 if(UpLo==(Upper|Lower))
318 count[StorageOrderMatch ? jp : ip]++;
319 else if(r==c)
320 count[ip]++;
321 else if(( UpLo==Lower && r>c) || ( UpLo==Upper && r<c))
322 {
323 count[ip]++;
324 count[jp]++;
325 }
326 }
327 }
328 Index nnz = count.sum();
329
330 // reserve space
331 dest.resizeNonZeros(nnz);
332 dest.outerIndexPtr()[0] = 0;
333 for(Index j=0; j<size; ++j)
334 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
335 for(Index j=0; j<size; ++j)
336 count[j] = dest.outerIndexPtr()[j];
337
338 // copy data
339 for(Index j = 0; j<size; ++j)
340 {
341 for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
342 {
343 Index i = it.index();
344 Index r = it.row();
345 Index c = it.col();
346
347 Index jp = perm ? perm[j] : j;
348 Index ip = perm ? perm[i] : i;
349
350 if(UpLo==(Upper|Lower))
351 {
352 Index k = count[StorageOrderMatch ? jp : ip]++;
353 dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp;
354 dest.valuePtr()[k] = it.value();
355 }
356 else if(r==c)
357 {
358 Index k = count[ip]++;
359 dest.innerIndexPtr()[k] = ip;
360 dest.valuePtr()[k] = it.value();
361 }
362 else if(( (UpLo&Lower)==Lower && r>c) || ( (UpLo&Upper)==Upper && r<c))
363 {
364 if(!StorageOrderMatch)
365 std::swap(ip,jp);
366 Index k = count[jp]++;
367 dest.innerIndexPtr()[k] = ip;
368 dest.valuePtr()[k] = it.value();
369 k = count[ip]++;
370 dest.innerIndexPtr()[k] = jp;
371 dest.valuePtr()[k] = internal::conj(it.value());
372 }
373 }
374 }
375}
376
377template<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder>
378void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
379{
380 typedef typename MatrixType::Index Index;
381 typedef typename MatrixType::Scalar Scalar;
382 SparseMatrix<Scalar,DstOrder,Index>& dest(_dest.derived());
383 typedef Matrix<Index,Dynamic,1> VectorI;
384 enum {
385 SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor,
386 StorageOrderMatch = int(SrcOrder) == int(DstOrder),
387 DstUpLo = DstOrder==RowMajor ? (_DstUpLo==Upper ? Lower : Upper) : _DstUpLo,
388 SrcUpLo = SrcOrder==RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo
389 };
390
391 Index size = mat.rows();
392 VectorI count(size);
393 count.setZero();
394 dest.resize(size,size);
395 for(Index j = 0; j<size; ++j)
396 {
397 Index jp = perm ? perm[j] : j;
398 for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
399 {
400 Index i = it.index();
401 if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j))
402 continue;
403
404 Index ip = perm ? perm[i] : i;
405 count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
406 }
407 }
408 dest.outerIndexPtr()[0] = 0;
409 for(Index j=0; j<size; ++j)
410 dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
411 dest.resizeNonZeros(dest.outerIndexPtr()[size]);
412 for(Index j=0; j<size; ++j)
413 count[j] = dest.outerIndexPtr()[j];
414
415 for(Index j = 0; j<size; ++j)
416 {
417
418 for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
419 {
420 Index i = it.index();
421 if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j))
422 continue;
423
424 Index jp = perm ? perm[j] : j;
425 Index ip = perm? perm[i] : i;
426
427 Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
428 dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp);
429
430 if(!StorageOrderMatch) std::swap(ip,jp);
431 if( ((int(DstUpLo)==int(Lower) && ip<jp) || (int(DstUpLo)==int(Upper) && ip>jp)))
432 dest.valuePtr()[k] = conj(it.value());
433 else
434 dest.valuePtr()[k] = it.value();
435 }
436 }
437}
438
439}
440
441template<typename MatrixType,int UpLo>
442class SparseSymmetricPermutationProduct
443 : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
444{
445 public:
446 typedef typename MatrixType::Scalar Scalar;
447 typedef typename MatrixType::Index Index;
448 protected:
449 typedef PermutationMatrix<Dynamic,Dynamic,Index> Perm;
450 public:
451 typedef Matrix<Index,Dynamic,1> VectorI;
452 typedef typename MatrixType::Nested MatrixTypeNested;
453 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
454
455 SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm)
456 : m_matrix(mat), m_perm(perm)
457 {}
458
459 inline Index rows() const { return m_matrix.rows(); }
460 inline Index cols() const { return m_matrix.cols(); }
461
462 template<typename DestScalar, int Options, typename DstIndex>
463 void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest) const
464 {
465 internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
466 }
467
468 template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const
469 {
470 internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
471 }
472
473 protected:
474 MatrixTypeNested m_matrix;
475 const Perm& m_perm;
476
477};
478
479} // end namespace Eigen
480
481#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H
Eigen::DenseBase::Index
internal::traits< Derived >::Index Index
The type of indices.
Definition DenseBase.h:51
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:50
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:129
Eigen::PermutationMatrix
Permutation matrix.
Definition PermutationMatrix.h:284
Eigen::SparseMatrixBase
Base class of any sparse matrices or sparse expressions.
Definition SparseMatrixBase.h:27
Eigen::SparseMatrix
A versatible sparse matrix representation.
Definition SparseMatrix.h:87
Eigen::SparseSelfAdjointView
Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
Definition SparseSelfAdjointView.h:51
Eigen::SparseSelfAdjointView::rankUpdate
SparseSelfAdjointView & rankUpdate(const SparseMatrixBase< DerivedU > &u, Scalar alpha=Scalar(1))
Eigen::SparseSelfAdjointView::operator*
SparseSelfAdjointTimeDenseProduct< MatrixType, OtherDerived, UpLo > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition SparseSelfAdjointView.h:75
Eigen::SparseSelfAdjointView::twistedBy
SparseSymmetricPermutationProduct< _MatrixTypeNested, UpLo > twistedBy(const PermutationMatrix< Dynamic, Dynamic, Index > &perm) const
Definition SparseSelfAdjointView.h:114
Eigen::SparseSelfAdjointView::operator*
friend DenseTimeSparseSelfAdjointProduct< OtherDerived, MatrixType, UpLo > operator*(const MatrixBase< OtherDerived > &lhs, const SparseSelfAdjointView &rhs)
Definition SparseSelfAdjointView.h:83
Eigen::RowMajor
@ RowMajor
Definition Constants.h:259
Eigen::ColMajor
@ ColMajor
Definition Constants.h:257
Eigen::Upper
@ Upper
Definition Constants.h:164
Eigen::Lower
@ Lower
Definition Constants.h:162
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition Constants.h:48
Eigen::internal
Definition LDLT.h:18
Eigen::EigenBase
Definition EigenBase.h:27