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