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
13namespace Eigen {
14
15namespace internal {
16template<typename MatrixType, int Size>
17struct 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
35template<typename MatrixType, int Size>
36class 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
107template<typename _Scalar, int _Options, typename _Index, int Size>
108class 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
288template<typename Derived>
289SparseInnerVectorSet<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
297template<typename Derived>
298const 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
305template<typename Derived>
306SparseInnerVectorSet<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
314template<typename Derived>
315const 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
324template<typename Derived>
325SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(Index outer)
326{ return SparseInnerVectorSet<Derived,1>(derived(), outer); }
327
331template<typename Derived>
332const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(Index outer) const
333{ return SparseInnerVectorSet<Derived,1>(derived(), outer); }
334
336template<typename Derived>
337SparseInnerVectorSet<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
345template<typename Derived>
346const 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
353template<typename Derived>
354SparseInnerVectorSet<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
362template<typename Derived>
363const 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
374template<typename Derived>
375SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize)
376{ return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
377
381template<typename Derived>
382const 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
Eigen::SparseMatrixBase::innerVectors
SparseInnerVectorSet< Derived, Dynamic > innerVectors(Index outerStart, Index outerSize)
Definition SparseBlock.h:375
Eigen::SparseMatrixBase::col
SparseInnerVectorSet< Derived, 1 > col(Index j)
Definition SparseBlock.h:306
Eigen::SparseMatrixBase::innerVector
SparseInnerVectorSet< Derived, 1 > innerVector(Index outer)
Definition SparseBlock.h:325
Eigen::SparseMatrixBase::middleCols
SparseInnerVectorSet< Derived, Dynamic > middleCols(Index start, Index size)
Definition SparseBlock.h:354
Eigen::SparseMatrixBase::middleRows
SparseInnerVectorSet< Derived, Dynamic > middleRows(Index start, Index size)
Definition SparseBlock.h:337
Eigen::SparseMatrixBase::row
SparseInnerVectorSet< Derived, 1 > row(Index i)
Definition SparseBlock.h:289
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition Constants.h:48
Eigen::internal
Definition LDLT.h:18
Eigen::EigenBase::size
Index size() const
Definition EigenBase.h:49
Eigen::EigenBase::derived
Derived & derived()
Definition EigenBase.h:34