SparseSelfAdjointView.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_SPARSE_SELFADJOINTVIEW_H
26 #define EIGEN_SPARSE_SELFADJOINTVIEW_H
27 
28 namespace Eigen {
29 
44 template<typename Lhs, typename Rhs, int UpLo>
45 class SparseSelfAdjointTimeDenseProduct;
46 
47 template<typename Lhs, typename Rhs, int UpLo>
48 class DenseTimeSparseSelfAdjointProduct;
49 
50 namespace internal {
51 
52 template<typename MatrixType, unsigned int UpLo>
53 struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> {
54 };
55 
56 template<int SrcUpLo,int DstUpLo,typename MatrixType,int DestOrder>
57 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
58 
59 template<int UpLo,typename MatrixType,int DestOrder>
60 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm = 0);
61 
62 }
63 
64 template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
65  : public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> >
66 {
67  public:
68 
69  typedef typename MatrixType::Scalar Scalar;
70  typedef typename MatrixType::Index Index;
72  typedef typename MatrixType::Nested MatrixTypeNested;
73  typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
74 
75  inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix)
76  {
77  eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices");
78  }
79 
80  inline Index rows() const { return m_matrix.rows(); }
81  inline Index cols() const { return m_matrix.cols(); }
82 
84  const _MatrixTypeNested& matrix() const { return m_matrix; }
85  _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); }
86 
88  template<typename OtherDerived>
91  {
93  }
94 
96  template<typename OtherDerived> friend
99  {
101  }
102 
111  template<typename DerivedU>
113 
115  template<typename DestScalar,int StorageOrder> void evalTo(SparseMatrix<DestScalar,StorageOrder,Index>& _dest) const
116  {
117  internal::permute_symm_to_fullsymm<UpLo>(m_matrix, _dest);
118  }
119 
120  template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar,ColMajor,Index>& _dest) const
121  {
122  // TODO directly evaluate into _dest;
124  internal::permute_symm_to_fullsymm<UpLo>(m_matrix, tmp);
125  _dest = tmp;
126  }
127 
130  {
132  }
133 
134  template<typename SrcMatrixType,int SrcUpLo>
136  {
137  permutedMatrix.evalTo(*this);
138  return *this;
139  }
140 
141 
143  {
145  return *this = src.twistedBy(pnull);
146  }
147 
148  template<typename SrcMatrixType,unsigned int SrcUpLo>
150  {
152  return *this = src.twistedBy(pnull);
153  }
154 
155 
156  // const SparseLLT<PlainObject, UpLo> llt() const;
157  // const SparseLDLT<PlainObject, UpLo> ldlt() const;
158 
159  protected:
160 
161  typename MatrixType::Nested m_matrix;
164 };
165 
166 /***************************************************************************
167 * Implementation of SparseMatrixBase methods
168 ***************************************************************************/
169 
170 template<typename Derived>
171 template<unsigned int UpLo>
173 {
174  return derived();
175 }
176 
177 template<typename Derived>
178 template<unsigned int UpLo>
180 {
181  return derived();
182 }
183 
184 /***************************************************************************
185 * Implementation of SparseSelfAdjointView methods
186 ***************************************************************************/
187 
188 template<typename MatrixType, unsigned int UpLo>
189 template<typename DerivedU>
192 {
194  if(alpha==Scalar(0))
195  m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
196  else
197  m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>();
198 
199  return *this;
200 }
201 
202 /***************************************************************************
203 * Implementation of sparse self-adjoint time dense matrix
204 ***************************************************************************/
205 
206 namespace internal {
207 template<typename Lhs, typename Rhs, int UpLo>
208 struct traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> >
209  : traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
210 {
211  typedef Dense StorageKind;
212 };
213 }
214 
215 template<typename Lhs, typename Rhs, int UpLo>
217  : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
218 {
219  public:
221 
222  SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
223  {}
224 
225  template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
226  {
227  // TODO use alpha
228  eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry");
229  typedef typename internal::remove_all<Lhs>::type _Lhs;
230  typedef typename internal::remove_all<Rhs>::type _Rhs;
231  typedef typename _Lhs::InnerIterator LhsInnerIterator;
232  enum {
233  LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
234  ProcessFirstHalf =
235  ((UpLo&(Upper|Lower))==(Upper|Lower))
236  || ( (UpLo&Upper) && !LhsIsRowMajor)
237  || ( (UpLo&Lower) && LhsIsRowMajor),
238  ProcessSecondHalf = !ProcessFirstHalf
239  };
240  for (Index j=0; j<m_lhs.outerSize(); ++j)
241  {
242  LhsInnerIterator i(m_lhs,j);
243  if (ProcessSecondHalf)
244  {
245  while (i && i.index()<j) ++i;
246  if(i && i.index()==j)
247  {
248  dest.row(j) += i.value() * m_rhs.row(j);
249  ++i;
250  }
251  }
252  for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
253  {
254  Index a = LhsIsRowMajor ? j : i.index();
255  Index b = LhsIsRowMajor ? i.index() : j;
256  typename Lhs::Scalar v = i.value();
257  dest.row(a) += (v) * m_rhs.row(b);
258  dest.row(b) += internal::conj(v) * m_rhs.row(a);
259  }
260  if (ProcessFirstHalf && i && (i.index()==j))
261  dest.row(j) += i.value() * m_rhs.row(j);
262  }
263  }
264 
265  private:
267 };
268 
269 namespace internal {
270 template<typename Lhs, typename Rhs, int UpLo>
271 struct traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
272  : traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
273 {};
274 }
275 
276 template<typename Lhs, typename Rhs, int UpLo>
278  : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
279 {
280  public:
282 
283  DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
284  {}
285 
286  template<typename Dest> void scaleAndAddTo(Dest& /*dest*/, Scalar /*alpha*/) const
287  {
288  // TODO
289  }
290 
291  private:
293 };
294 
295 /***************************************************************************
296 * Implementation of symmetric copies and permutations
297 ***************************************************************************/
298 namespace internal {
299 
300 template<typename MatrixType, int UpLo>
301 struct traits<SparseSymmetricPermutationProduct<MatrixType,UpLo> > : traits<MatrixType> {
302 };
303 
304 template<int UpLo,typename MatrixType,int DestOrder>
305 void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DestOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
306 {
307  typedef typename MatrixType::Index Index;
308  typedef typename MatrixType::Scalar Scalar;
310  typedef Matrix<Index,Dynamic,1> VectorI;
311 
312  Dest& dest(_dest.derived());
313  enum {
314  StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
315  };
316 
317  Index size = mat.rows();
318  VectorI count;
319  count.resize(size);
320  count.setZero();
321  dest.resize(size,size);
322  for(Index j = 0; j<size; ++j)
323  {
324  Index jp = perm ? perm[j] : j;
325  for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
326  {
327  Index i = it.index();
328  Index r = it.row();
329  Index c = it.col();
330  Index ip = perm ? perm[i] : i;
331  if(UpLo==(Upper|Lower))
332  count[StorageOrderMatch ? jp : ip]++;
333  else if(r==c)
334  count[ip]++;
335  else if(( UpLo==Lower && r>c) || ( UpLo==Upper && r<c))
336  {
337  count[ip]++;
338  count[jp]++;
339  }
340  }
341  }
342  Index nnz = count.sum();
343 
344  // reserve space
345  dest.resizeNonZeros(nnz);
346  dest.outerIndexPtr()[0] = 0;
347  for(Index j=0; j<size; ++j)
348  dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
349  for(Index j=0; j<size; ++j)
350  count[j] = dest.outerIndexPtr()[j];
351 
352  // copy data
353  for(Index j = 0; j<size; ++j)
354  {
355  for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
356  {
357  Index i = it.index();
358  Index r = it.row();
359  Index c = it.col();
360 
361  Index jp = perm ? perm[j] : j;
362  Index ip = perm ? perm[i] : i;
363 
364  if(UpLo==(Upper|Lower))
365  {
366  Index k = count[StorageOrderMatch ? jp : ip]++;
367  dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp;
368  dest.valuePtr()[k] = it.value();
369  }
370  else if(r==c)
371  {
372  Index k = count[ip]++;
373  dest.innerIndexPtr()[k] = ip;
374  dest.valuePtr()[k] = it.value();
375  }
376  else if(( (UpLo&Lower)==Lower && r>c) || ( (UpLo&Upper)==Upper && r<c))
377  {
378  if(!StorageOrderMatch)
379  std::swap(ip,jp);
380  Index k = count[jp]++;
381  dest.innerIndexPtr()[k] = ip;
382  dest.valuePtr()[k] = it.value();
383  k = count[ip]++;
384  dest.innerIndexPtr()[k] = jp;
385  dest.valuePtr()[k] = internal::conj(it.value());
386  }
387  }
388  }
389 }
390 
391 template<int _SrcUpLo,int _DstUpLo,typename MatrixType,int DstOrder>
392 void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixType::Scalar,DstOrder,typename MatrixType::Index>& _dest, const typename MatrixType::Index* perm)
393 {
394  typedef typename MatrixType::Index Index;
395  typedef typename MatrixType::Scalar Scalar;
397  typedef Matrix<Index,Dynamic,1> VectorI;
398  enum {
399  SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor,
400  StorageOrderMatch = int(SrcOrder) == int(DstOrder),
401  DstUpLo = DstOrder==RowMajor ? (_DstUpLo==Upper ? Lower : Upper) : _DstUpLo,
402  SrcUpLo = SrcOrder==RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo
403  };
404 
405  Index size = mat.rows();
406  VectorI count(size);
407  count.setZero();
408  dest.resize(size,size);
409  for(Index j = 0; j<size; ++j)
410  {
411  Index jp = perm ? perm[j] : j;
412  for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
413  {
414  Index i = it.index();
415  if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j))
416  continue;
417 
418  Index ip = perm ? perm[i] : i;
419  count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
420  }
421  }
422  dest.outerIndexPtr()[0] = 0;
423  for(Index j=0; j<size; ++j)
424  dest.outerIndexPtr()[j+1] = dest.outerIndexPtr()[j] + count[j];
425  dest.resizeNonZeros(dest.outerIndexPtr()[size]);
426  for(Index j=0; j<size; ++j)
427  count[j] = dest.outerIndexPtr()[j];
428 
429  for(Index j = 0; j<size; ++j)
430  {
431 
432  for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
433  {
434  Index i = it.index();
435  if((int(SrcUpLo)==int(Lower) && i<j) || (int(SrcUpLo)==int(Upper) && i>j))
436  continue;
437 
438  Index jp = perm ? perm[j] : j;
439  Index ip = perm? perm[i] : i;
440 
441  Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
442  dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp);
443 
444  if(!StorageOrderMatch) std::swap(ip,jp);
445  if( ((int(DstUpLo)==int(Lower) && ip<jp) || (int(DstUpLo)==int(Upper) && ip>jp)))
446  dest.valuePtr()[k] = conj(it.value());
447  else
448  dest.valuePtr()[k] = it.value();
449  }
450  }
451 }
452 
453 }
454 
455 template<typename MatrixType,int UpLo>
457  : public EigenBase<SparseSymmetricPermutationProduct<MatrixType,UpLo> >
458 {
459  public:
460  typedef typename MatrixType::Scalar Scalar;
461  typedef typename MatrixType::Index Index;
462  protected:
464  public:
466  typedef typename MatrixType::Nested MatrixTypeNested;
467  typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
468 
469  SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm)
470  : m_matrix(mat), m_perm(perm)
471  {}
472 
473  inline Index rows() const { return m_matrix.rows(); }
474  inline Index cols() const { return m_matrix.cols(); }
475 
476  template<typename DestScalar, int Options, typename DstIndex>
478  {
479  internal::permute_symm_to_fullsymm<UpLo>(m_matrix,_dest,m_perm.indices().data());
480  }
481 
482  template<typename DestType,unsigned int DestUpLo> void evalTo(SparseSelfAdjointView<DestType,DestUpLo>& dest) const
483  {
484  internal::permute_symm_to_symm<UpLo,DestUpLo>(m_matrix,dest.matrix(),m_perm.indices().data());
485  }
486 
487  protected:
489  const Perm& m_perm;
490 
491 };
492 
493 } // end namespace Eigen
494 
495 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H