SparseMatrixBase.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2011 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_SPARSEMATRIXBASE_H
11 #define EIGEN_SPARSEMATRIXBASE_H
12 
13 namespace Eigen {
14 
26 template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
27 {
28  public:
29 
30  typedef typename internal::traits<Derived>::Scalar Scalar;
31  typedef typename internal::packet_traits<Scalar>::type PacketScalar;
32  typedef typename internal::traits<Derived>::StorageKind StorageKind;
33  typedef typename internal::traits<Derived>::Index Index;
34  typedef typename internal::add_const_on_value_type_if_arithmetic<
35  typename internal::packet_traits<Scalar>::type
36  >::type PacketReturnType;
37 
39  typedef EigenBase<Derived> Base;
40 
41  template<typename OtherDerived>
42  Derived& operator=(const EigenBase<OtherDerived> &other)
43  {
44  other.derived().evalTo(derived());
45  return derived();
46  }
47 
48  enum {
49 
50  RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
56  ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
63  SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime,
64  internal::traits<Derived>::ColsAtCompileTime>::ret),
69  MaxRowsAtCompileTime = RowsAtCompileTime,
70  MaxColsAtCompileTime = ColsAtCompileTime,
71 
72  MaxSizeAtCompileTime = (internal::size_at_compile_time<MaxRowsAtCompileTime,
73  MaxColsAtCompileTime>::ret),
74 
81  Flags = internal::traits<Derived>::Flags,
86  CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
91  IsRowMajor = Flags&RowMajorBit ? 1 : 0,
92 
93  #ifndef EIGEN_PARSED_BY_DOXYGEN
94  _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC
95  #endif
96  };
97 
99  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
102  >::type AdjointReturnType;
103 
104 
106 
107 
108 #ifndef EIGEN_PARSED_BY_DOXYGEN
115  typedef typename NumTraits<Scalar>::Real RealScalar;
116 
119  typedef typename internal::conditional<_HasDirectAccess, const Scalar&, Scalar>::type CoeffReturnType;
120 
123 
125  typedef Matrix<Scalar,EIGEN_SIZE_MAX(RowsAtCompileTime,ColsAtCompileTime),
126  EIGEN_SIZE_MAX(RowsAtCompileTime,ColsAtCompileTime)> SquareMatrixType;
127 
128  inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
129  inline Derived& derived() { return *static_cast<Derived*>(this); }
130  inline Derived& const_cast_derived() const
131  { return *static_cast<Derived*>(const_cast<SparseMatrixBase*>(this)); }
132 #endif // not EIGEN_PARSED_BY_DOXYGEN
133 
134 #define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::SparseMatrixBase
135 # include "../plugins/CommonCwiseUnaryOps.h"
136 # include "../plugins/CommonCwiseBinaryOps.h"
137 # include "../plugins/MatrixCwiseUnaryOps.h"
138 # include "../plugins/MatrixCwiseBinaryOps.h"
139 # ifdef EIGEN_SPARSEMATRIXBASE_PLUGIN
140 # include EIGEN_SPARSEMATRIXBASE_PLUGIN
141 # endif
142 # undef EIGEN_CURRENT_STORAGE_BASE_CLASS
143 #undef EIGEN_CURRENT_STORAGE_BASE_CLASS
144 
145 
147  inline Index rows() const { return derived().rows(); }
149  inline Index cols() const { return derived().cols(); }
152  inline Index size() const { return rows() * cols(); }
155  inline Index nonZeros() const { return derived().nonZeros(); }
160  inline bool isVector() const { return rows()==1 || cols()==1; }
163  Index outerSize() const { return (int(Flags)&RowMajorBit) ? this->rows() : this->cols(); }
166  Index innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); }
168  bool isRValue() const { return m_isRValue; }
169  Derived& markAsRValue() { m_isRValue = true; return derived(); }
170 
171  SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ }
172 
173 
174  template<typename OtherDerived>
175  Derived& operator=(const ReturnByValue<OtherDerived>& other)
176  {
177  other.evalTo(derived());
178  return derived();
179  }
180 
181 
182  template<typename OtherDerived>
183  inline Derived& operator=(const SparseMatrixBase<OtherDerived>& other)
184  {
185  return assign(other.derived());
186  }
187 
188  inline Derived& operator=(const Derived& other)
189  {
190 // if (other.isRValue())
191 // derived().swap(other.const_cast_derived());
192 // else
193  return assign(other.derived());
194  }
195 
196  protected:
197 
198  template<typename OtherDerived>
199  inline Derived& assign(const OtherDerived& other)
200  {
201  const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
202  const Index outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols();
203  if ((!transpose) && other.isRValue())
204  {
205  // eval without temporary
206  derived().resize(other.rows(), other.cols());
207  derived().setZero();
208  derived().reserve((std::max)(this->rows(),this->cols())*2);
209  for (Index j=0; j<outerSize; ++j)
210  {
211  derived().startVec(j);
212  for (typename OtherDerived::InnerIterator it(other, j); it; ++it)
213  {
214  Scalar v = it.value();
215  derived().insertBackByOuterInner(j,it.index()) = v;
216  }
217  }
218  derived().finalize();
219  }
220  else
221  {
222  assignGeneric(other);
223  }
224  return derived();
225  }
226 
227  template<typename OtherDerived>
228  inline void assignGeneric(const OtherDerived& other)
229  {
230  //const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
231  eigen_assert(( ((internal::traits<Derived>::SupportedAccessPatterns&OuterRandomAccessPattern)==OuterRandomAccessPattern) ||
232  (!((Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit)))) &&
233  "the transpose operation is supposed to be handled in SparseMatrix::operator=");
234 
235  enum { Flip = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit) };
236 
237  const Index outerSize = other.outerSize();
238  //typedef typename internal::conditional<transpose, LinkedVectorMatrix<Scalar,Flags&RowMajorBit>, Derived>::type TempType;
239  // thanks to shallow copies, we always eval to a tempary
240  Derived temp(other.rows(), other.cols());
241 
242  temp.reserve((std::max)(this->rows(),this->cols())*2);
243  for (Index j=0; j<outerSize; ++j)
244  {
245  temp.startVec(j);
246  for (typename OtherDerived::InnerIterator it(other.derived(), j); it; ++it)
247  {
248  Scalar v = it.value();
249  temp.insertBackByOuterInner(Flip?it.index():j,Flip?j:it.index()) = v;
250  }
251  }
252  temp.finalize();
253 
254  derived() = temp.markAsRValue();
255  }
256 
257  public:
258 
259  template<typename Lhs, typename Rhs>
260  inline Derived& operator=(const SparseSparseProduct<Lhs,Rhs>& product);
261 
262  friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m)
263  {
264  typedef typename Derived::Nested Nested;
265  typedef typename internal::remove_all<Nested>::type NestedCleaned;
266 
267  if (Flags&RowMajorBit)
268  {
269  const Nested nm(m.derived());
270  for (Index row=0; row<nm.outerSize(); ++row)
271  {
272  Index col = 0;
273  for (typename NestedCleaned::InnerIterator it(nm.derived(), row); it; ++it)
274  {
275  for ( ; col<it.index(); ++col)
276  s << "0 ";
277  s << it.value() << " ";
278  ++col;
279  }
280  for ( ; col<m.cols(); ++col)
281  s << "0 ";
282  s << std::endl;
283  }
284  }
285  else
286  {
287  const Nested nm(m.derived());
288  if (m.cols() == 1) {
289  Index row = 0;
290  for (typename NestedCleaned::InnerIterator it(nm.derived(), 0); it; ++it)
291  {
292  for ( ; row<it.index(); ++row)
293  s << "0" << std::endl;
294  s << it.value() << std::endl;
295  ++row;
296  }
297  for ( ; row<m.rows(); ++row)
298  s << "0" << std::endl;
299  }
300  else
301  {
302  SparseMatrix<Scalar, RowMajorBit> trans = m;
303  s << static_cast<const SparseMatrixBase<SparseMatrix<Scalar, RowMajorBit> >&>(trans);
304  }
305  }
306  return s;
307  }
308 
309  template<typename OtherDerived>
310  Derived& operator+=(const SparseMatrixBase<OtherDerived>& other);
311  template<typename OtherDerived>
312  Derived& operator-=(const SparseMatrixBase<OtherDerived>& other);
313 
314  Derived& operator*=(const Scalar& other);
315  Derived& operator/=(const Scalar& other);
316 
317  #define EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE \
318  CwiseBinaryOp< \
319  internal::scalar_product_op< \
320  typename internal::scalar_product_traits< \
321  typename internal::traits<Derived>::Scalar, \
322  typename internal::traits<OtherDerived>::Scalar \
323  >::ReturnType \
324  >, \
325  Derived, \
326  OtherDerived \
327  >
328 
329  template<typename OtherDerived>
330  EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
331  cwiseProduct(const MatrixBase<OtherDerived> &other) const;
332 
333  // sparse * sparse
334  template<typename OtherDerived>
335  const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
336  operator*(const SparseMatrixBase<OtherDerived> &other) const;
337 
338  // sparse * diagonal
339  template<typename OtherDerived>
340  const SparseDiagonalProduct<Derived,OtherDerived>
341  operator*(const DiagonalBase<OtherDerived> &other) const;
342 
343  // diagonal * sparse
344  template<typename OtherDerived> friend
345  const SparseDiagonalProduct<OtherDerived,Derived>
346  operator*(const DiagonalBase<OtherDerived> &lhs, const SparseMatrixBase& rhs)
347  { return SparseDiagonalProduct<OtherDerived,Derived>(lhs.derived(), rhs.derived()); }
348 
350  template<typename OtherDerived> friend
351  const typename DenseSparseProductReturnType<OtherDerived,Derived>::Type
352  operator*(const MatrixBase<OtherDerived>& lhs, const Derived& rhs)
353  { return typename DenseSparseProductReturnType<OtherDerived,Derived>::Type(lhs.derived(),rhs); }
354 
356  template<typename OtherDerived>
357  const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
358  operator*(const MatrixBase<OtherDerived> &other) const;
359 
361  SparseSymmetricPermutationProduct<Derived,Upper|Lower> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const
362  {
363  return SparseSymmetricPermutationProduct<Derived,Upper|Lower>(derived(), perm);
364  }
365 
366  template<typename OtherDerived>
367  Derived& operator*=(const SparseMatrixBase<OtherDerived>& other);
368 
369  #ifdef EIGEN2_SUPPORT
370  // deprecated
371  template<typename OtherDerived>
373  solveTriangular(const MatrixBase<OtherDerived>& other) const;
374 
375  // deprecated
376  template<typename OtherDerived>
377  void solveTriangularInPlace(MatrixBase<OtherDerived>& other) const;
378  #endif // EIGEN2_SUPPORT
379 
380  template<int Mode>
381  inline const SparseTriangularView<Derived, Mode> triangularView() const;
382 
383  template<unsigned int UpLo> inline const SparseSelfAdjointView<Derived, UpLo> selfadjointView() const;
384  template<unsigned int UpLo> inline SparseSelfAdjointView<Derived, UpLo> selfadjointView();
385 
386  template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const;
387  template<typename OtherDerived> Scalar dot(const SparseMatrixBase<OtherDerived>& other) const;
388  RealScalar squaredNorm() const;
389  RealScalar norm() const;
390 
391  Transpose<Derived> transpose() { return derived(); }
392  const Transpose<const Derived> transpose() const { return derived(); }
393  const AdjointReturnType adjoint() const { return transpose(); }
394 
395  // sub-vector
396  SparseInnerVectorSet<Derived,1> row(Index i);
397  const SparseInnerVectorSet<Derived,1> row(Index i) const;
398  SparseInnerVectorSet<Derived,1> col(Index j);
399  const SparseInnerVectorSet<Derived,1> col(Index j) const;
400  SparseInnerVectorSet<Derived,1> innerVector(Index outer);
401  const SparseInnerVectorSet<Derived,1> innerVector(Index outer) const;
402 
403  // set of sub-vectors
404  SparseInnerVectorSet<Derived,Dynamic> subrows(Index start, Index size);
405  const SparseInnerVectorSet<Derived,Dynamic> subrows(Index start, Index size) const;
406  SparseInnerVectorSet<Derived,Dynamic> subcols(Index start, Index size);
407  const SparseInnerVectorSet<Derived,Dynamic> subcols(Index start, Index size) const;
408 
409  SparseInnerVectorSet<Derived,Dynamic> middleRows(Index start, Index size);
410  const SparseInnerVectorSet<Derived,Dynamic> middleRows(Index start, Index size) const;
411  SparseInnerVectorSet<Derived,Dynamic> middleCols(Index start, Index size);
412  const SparseInnerVectorSet<Derived,Dynamic> middleCols(Index start, Index size) const;
413  SparseInnerVectorSet<Derived,Dynamic> innerVectors(Index outerStart, Index outerSize);
414  const SparseInnerVectorSet<Derived,Dynamic> innerVectors(Index outerStart, Index outerSize) const;
415 
417  template<typename DenseDerived>
418  void evalTo(MatrixBase<DenseDerived>& dst) const
419  {
420  dst.setZero();
421  for (Index j=0; j<outerSize(); ++j)
422  for (typename Derived::InnerIterator i(derived(),j); i; ++i)
423  dst.coeffRef(i.row(),i.col()) = i.value();
424  }
425 
426  Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> toDense() const
427  {
428  return derived();
429  }
430 
431  template<typename OtherDerived>
432  bool isApprox(const SparseMatrixBase<OtherDerived>& other,
433  RealScalar prec = NumTraits<Scalar>::dummy_precision()) const
434  { return toDense().isApprox(other.toDense(),prec); }
435 
436  template<typename OtherDerived>
437  bool isApprox(const MatrixBase<OtherDerived>& other,
438  RealScalar prec = NumTraits<Scalar>::dummy_precision()) const
439  { return toDense().isApprox(other,prec); }
440 
446  inline const typename internal::eval<Derived>::type eval() const
447  { return typename internal::eval<Derived>::type(derived()); }
448 
449  Scalar sum() const;
450 
451  protected:
452 
453  bool m_isRValue;
454 };
455 
456 } // end namespace Eigen
457 
458 #endif // EIGEN_SPARSEMATRIXBASE_H