SparseBlock.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-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_BLOCK_H
11 #define EIGEN_SPARSE_BLOCK_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 template<typename MatrixType, int Size>
17 struct traits<SparseInnerVectorSet<MatrixType, Size> >
18 {
19  typedef typename traits<MatrixType>::Scalar Scalar;
20  typedef typename traits<MatrixType>::Index Index;
21  typedef typename traits<MatrixType>::StorageKind StorageKind;
22  typedef MatrixXpr XprKind;
23  enum {
24  IsRowMajor = (int(MatrixType::Flags)&RowMajorBit)==RowMajorBit,
25  Flags = MatrixType::Flags,
26  RowsAtCompileTime = IsRowMajor ? Size : MatrixType::RowsAtCompileTime,
27  ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : Size,
28  MaxRowsAtCompileTime = RowsAtCompileTime,
29  MaxColsAtCompileTime = ColsAtCompileTime,
30  CoeffReadCost = MatrixType::CoeffReadCost
31  };
32 };
33 } // end namespace internal
34 
35 template<typename MatrixType, int Size>
36 class SparseInnerVectorSet : internal::no_assignment_operator,
37  public SparseMatrixBase<SparseInnerVectorSet<MatrixType, Size> >
38 {
39  public:
40 
41  enum { IsRowMajor = internal::traits<SparseInnerVectorSet>::IsRowMajor };
42 
43  EIGEN_SPARSE_PUBLIC_INTERFACE(SparseInnerVectorSet)
44  class InnerIterator: public MatrixType::InnerIterator
45  {
46  public:
47  inline InnerIterator(const SparseInnerVectorSet& xpr, Index outer)
48  : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
49  {}
50  inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
51  inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
52  protected:
53  Index m_outer;
54  };
55  class ReverseInnerIterator: public MatrixType::ReverseInnerIterator
56  {
57  public:
58  inline ReverseInnerIterator(const SparseInnerVectorSet& xpr, Index outer)
59  : MatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
60  {}
61  inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
62  inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
63  protected:
64  Index m_outer;
65  };
66 
67  inline SparseInnerVectorSet(const MatrixType& matrix, Index outerStart, Index outerSize)
68  : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
69  {
70  eigen_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
71  }
72 
73  inline SparseInnerVectorSet(const MatrixType& matrix, Index outer)
74  : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
75  {
76  eigen_assert(Size!=Dynamic);
77  eigen_assert( (outer>=0) && (outer<matrix.outerSize()) );
78  }
79 
80 // template<typename OtherDerived>
81 // inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
82 // {
83 // return *this;
84 // }
85 
86 // template<typename Sparse>
87 // inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
88 // {
89 // return *this;
90 // }
91 
92  EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
93  EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
94 
95  protected:
96 
97  const typename MatrixType::Nested m_matrix;
98  Index m_outerStart;
99  const internal::variable_if_dynamic<Index, Size> m_outerSize;
100 };
101 
102 
103 /***************************************************************************
104 * specialisation for SparseMatrix
105 ***************************************************************************/
106 
107 template<typename _Scalar, int _Options, typename _Index, int Size>
108 class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options, _Index>, Size>
109  : public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options, _Index>, Size> >
110 {
111  typedef SparseMatrix<_Scalar, _Options, _Index> MatrixType;
112  public:
113 
114  enum { IsRowMajor = internal::traits<SparseInnerVectorSet>::IsRowMajor };
115 
116  EIGEN_SPARSE_PUBLIC_INTERFACE(SparseInnerVectorSet)
117  class InnerIterator: public MatrixType::InnerIterator
118  {
119  public:
120  inline InnerIterator(const SparseInnerVectorSet& xpr, Index outer)
121  : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
122  {}
123  inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
124  inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
125  protected:
126  Index m_outer;
127  };
128  class ReverseInnerIterator: public MatrixType::ReverseInnerIterator
129  {
130  public:
131  inline ReverseInnerIterator(const SparseInnerVectorSet& xpr, Index outer)
132  : MatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
133  {}
134  inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
135  inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
136  protected:
137  Index m_outer;
138  };
139 
140  inline SparseInnerVectorSet(const MatrixType& matrix, Index outerStart, Index outerSize)
141  : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
142  {
143  eigen_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
144  }
145 
146  inline SparseInnerVectorSet(const MatrixType& matrix, Index outer)
147  : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
148  {
149  eigen_assert(Size==1);
150  eigen_assert( (outer>=0) && (outer<matrix.outerSize()) );
151  }
152 
153  template<typename OtherDerived>
154  inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
155  {
156  typedef typename internal::remove_all<typename MatrixType::Nested>::type _NestedMatrixType;
157  _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);;
158  // This assignement is slow if this vector set is not empty
159  // and/or it is not at the end of the nonzeros of the underlying matrix.
160 
161  // 1 - eval to a temporary to avoid transposition and/or aliasing issues
162  SparseMatrix<Scalar, IsRowMajor ? RowMajor : ColMajor, Index> tmp(other);
163 
164  // 2 - let's check whether there is enough allocated memory
165  Index nnz = tmp.nonZeros();
166  Index nnz_previous = nonZeros();
167  Index free_size = Index(matrix.data().allocatedSize()) + nnz_previous;
168  Index nnz_head = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart];
169  Index tail = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()];
170  Index nnz_tail = matrix.nonZeros() - tail;
171 
172  if(nnz>free_size)
173  {
174  // realloc manually to reduce copies
175  typename MatrixType::Storage newdata(m_matrix.nonZeros() - nnz_previous + nnz);
176 
177  std::memcpy(&newdata.value(0), &m_matrix.data().value(0), nnz_head*sizeof(Scalar));
178  std::memcpy(&newdata.index(0), &m_matrix.data().index(0), nnz_head*sizeof(Index));
179 
180  std::memcpy(&newdata.value(nnz_head), &tmp.data().value(0), nnz*sizeof(Scalar));
181  std::memcpy(&newdata.index(nnz_head), &tmp.data().index(0), nnz*sizeof(Index));
182 
183  std::memcpy(&newdata.value(nnz_head+nnz), &matrix.data().value(tail), nnz_tail*sizeof(Scalar));
184  std::memcpy(&newdata.index(nnz_head+nnz), &matrix.data().index(tail), nnz_tail*sizeof(Index));
185 
186  matrix.data().swap(newdata);
187  }
188  else
189  {
190  // no need to realloc, simply copy the tail at its respective position and insert tmp
191  matrix.data().resize(nnz_head + nnz + nnz_tail);
192 
193  if(nnz<nnz_previous)
194  {
195  std::memcpy(&matrix.data().value(nnz_head+nnz), &matrix.data().value(tail), nnz_tail*sizeof(Scalar));
196  std::memcpy(&matrix.data().index(nnz_head+nnz), &matrix.data().index(tail), nnz_tail*sizeof(Index));
197  }
198  else
199  {
200  for(Index i=nnz_tail-1; i>=0; --i)
201  {
202  matrix.data().value(nnz_head+nnz+i) = matrix.data().value(tail+i);
203  matrix.data().index(nnz_head+nnz+i) = matrix.data().index(tail+i);
204  }
205  }
206 
207  std::memcpy(&matrix.data().value(nnz_head), &tmp.data().value(0), nnz*sizeof(Scalar));
208  std::memcpy(&matrix.data().index(nnz_head), &tmp.data().index(0), nnz*sizeof(Index));
209  }
210 
211  // update outer index pointers
212  Index p = nnz_head;
213  for(Index k=0; k<m_outerSize.value(); ++k)
214  {
215  matrix.outerIndexPtr()[m_outerStart+k] = p;
216  p += tmp.innerVector(k).nonZeros();
217  }
218  std::ptrdiff_t offset = nnz - nnz_previous;
219  for(Index k = m_outerStart + m_outerSize.value(); k<=matrix.outerSize(); ++k)
220  {
221  matrix.outerIndexPtr()[k] += offset;
222  }
223 
224  return *this;
225  }
226 
227  inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other)
228  {
229  return operator=<SparseInnerVectorSet>(other);
230  }
231 
232  inline const Scalar* valuePtr() const
233  { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
234  inline Scalar* valuePtr()
235  { return m_matrix.const_cast_derived().valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
236 
237  inline const Index* innerIndexPtr() const
238  { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
239  inline Index* innerIndexPtr()
240  { return m_matrix.const_cast_derived().innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
241 
242  inline const Index* outerIndexPtr() const
243  { return m_matrix.outerIndexPtr() + m_outerStart; }
244  inline Index* outerIndexPtr()
245  { return m_matrix.const_cast_derived().outerIndexPtr() + m_outerStart; }
246 
247  Index nonZeros() const
248  {
249  if(m_matrix.isCompressed())
250  return std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()])
251  - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]);
252  else if(m_outerSize.value()==0)
253  return 0;
254  else
255  return Map<const Matrix<Index,Size,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum();
256  }
257 
258  const Scalar& lastCoeff() const
259  {
260  EIGEN_STATIC_ASSERT_VECTOR_ONLY(SparseInnerVectorSet);
261  eigen_assert(nonZeros()>0);
262  if(m_matrix.isCompressed())
263  return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1];
264  else
265  return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1];
266  }
267 
268 // template<typename Sparse>
269 // inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
270 // {
271 // return *this;
272 // }
273 
274  EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
275  EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
276 
277  protected:
278 
279  typename MatrixType::Nested m_matrix;
280  Index m_outerStart;
281  const internal::variable_if_dynamic<Index, Size> m_outerSize;
282 
283 };
284 
285 //----------
286 
288 template<typename Derived>
289 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(Index i)
290 {
291  EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
292  return innerVector(i);
293 }
294 
297 template<typename Derived>
298 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(Index i) const
299 {
300  EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
301  return innerVector(i);
302 }
303 
305 template<typename Derived>
306 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(Index i)
307 {
308  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
309  return innerVector(i);
310 }
311 
314 template<typename Derived>
315 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(Index i) const
316 {
317  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
318  return innerVector(i);
319 }
320 
324 template<typename Derived>
325 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(Index outer)
326 { return SparseInnerVectorSet<Derived,1>(derived(), outer); }
327 
331 template<typename Derived>
332 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(Index outer) const
333 { return SparseInnerVectorSet<Derived,1>(derived(), outer); }
334 
336 template<typename Derived>
337 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleRows(Index start, Index size)
338 {
339  EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
340  return innerVectors(start, size);
341 }
342 
345 template<typename Derived>
346 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleRows(Index start, Index size) const
347 {
348  EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
349  return innerVectors(start, size);
350 }
351 
353 template<typename Derived>
354 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleCols(Index start, Index size)
355 {
356  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
357  return innerVectors(start, size);
358 }
359 
362 template<typename Derived>
363 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleCols(Index start, Index size) const
364 {
365  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
366  return innerVectors(start, size);
367 }
368 
369 
370 
374 template<typename Derived>
375 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize)
376 { return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
377 
381 template<typename Derived>
382 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize) const
383 { return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
384 
385 } // end namespace Eigen
386 
387 #endif // EIGEN_SPARSE_BLOCK_H