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  // 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 
254 namespace internal {
255 template<typename Lhs, typename Rhs, int UpLo>
256 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
257  : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
258 {};
259 }
260 
261 template<typename Lhs, typename Rhs, int UpLo>
262 class 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 ***************************************************************************/
283 namespace internal {
284 
285 template<typename MatrixType, int UpLo>
286 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
287 };
288 
289 template<int UpLo,typename MatrixType,int DestOrder>
290 void 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 
376 template<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder>
377 void 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 
440 template<typename MatrixType,int UpLo>
441 class 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