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 
13 namespace Eigen {
14 
29 template<typename Lhs, typename Rhs, int UpLo>
30 class SparseSelfAdjointTimeDenseProduct;
31 
32 template<typename Lhs, typename Rhs, int UpLo>
33 class DenseTimeSparseSelfAdjointProduct;
34 
35 namespace internal {
36 
37 template<typename MatrixType, unsigned int UpLo>
38 struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> {
39 };
40 
41 template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder>
42 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
43 
44 template<int UpLo,typename MatrixType,int DestOrder>
45 void 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 
49 template<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;
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>
76  {
77  return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived());
78  }
79 
81  template<typename OtherDerived> friend
82  DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>
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 
155 template<typename Derived>
156 template<unsigned int UpLo>
157 const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const
158 {
159  return derived();
160 }
161 
162 template<typename Derived>
163 template<unsigned int UpLo>
164 SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView()
165 {
166  return derived();
167 }
168 
169 /***************************************************************************
170 * Implementation of SparseSelfAdjointView methods
171 ***************************************************************************/
172 
173 template<typename MatrixType, unsigned int UpLo>
174 template<typename DerivedU>
175 SparseSelfAdjointView<MatrixType,UpLo>&
177 {
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 
191 namespace internal {
192 template<typename Lhs, typename Rhs, int UpLo>
193 struct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> >
194  : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
195 {
196  typedef Dense StorageKind;
197 };
198 }
199 
200 template<typename Lhs, typename Rhs, int UpLo>
201 class 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 
255 namespace internal {
256 template<typename Lhs, typename Rhs, int UpLo>
257 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
258  : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
259 {};
260 }
261 
262 template<typename Lhs, typename Rhs, int UpLo>
263 class 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 ***************************************************************************/
284 namespace internal {
285 
286 template<typename MatrixType, int UpLo>
287 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
288 };
289 
290 template<int UpLo,typename MatrixType,int DestOrder>
291 void 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 
377 template<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder>
378 void 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 
441 template<typename MatrixType,int UpLo>
442 class 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