1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2015 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_COMPRESSED_BASE_H
11 #define EIGEN_SPARSE_COMPRESSED_BASE_H
12 
13 namespace Eigen {
14 
15 template<typename Derived> class SparseCompressedBase;
16 
17 namespace internal {
18 
19 template<typename Derived>
20 struct traits<SparseCompressedBase<Derived> > : traits<Derived>
21 {};
22 
23 } // end namespace internal
24 
25 /** \ingroup SparseCore_Module
26   * \class SparseCompressedBase
27   * \brief Common base class for sparse [compressed]-{row|column}-storage format.
28   *
29   * This class defines the common interface for all derived classes implementing the compressed sparse storage format, such as:
30   *  - SparseMatrix
31   *  - Ref<SparseMatrixType,Options>
32   *  - Map<SparseMatrixType>
33   *
34   */
35 template<typename Derived>
36 class SparseCompressedBase
37   : public SparseMatrixBase<Derived>
38 {
39   public:
40     typedef SparseMatrixBase<Derived> Base;
41     EIGEN_SPARSE_PUBLIC_INTERFACE(SparseCompressedBase)
42     using Base::operator=;
43     using Base::IsRowMajor;
44 
45     class InnerIterator;
46     class ReverseInnerIterator;
47 
48   protected:
49     typedef typename Base::IndexVector IndexVector;
50     Eigen::Map<IndexVector> innerNonZeros() { return Eigen::Map<IndexVector>(innerNonZeroPtr(), isCompressed()?0:derived().outerSize()); }
51     const  Eigen::Map<const IndexVector> innerNonZeros() const { return Eigen::Map<const IndexVector>(innerNonZeroPtr(), isCompressed()?0:derived().outerSize()); }
52 
53   public:
54 
55     /** \returns the number of non zero coefficients */
56     inline Index nonZeros() const
57     {
58       if(Derived::IsVectorAtCompileTime && outerIndexPtr()==0)
59         return derived().nonZeros();
60       else if(isCompressed())
61         return outerIndexPtr()[derived().outerSize()]-outerIndexPtr()[0];
62       else if(derived().outerSize()==0)
63         return 0;
64       else
65         return innerNonZeros().sum();
66     }
67 
68     /** \returns a const pointer to the array of values.
69       * This function is aimed at interoperability with other libraries.
70       * \sa innerIndexPtr(), outerIndexPtr() */
71     inline const Scalar* valuePtr() const { return derived().valuePtr(); }
72     /** \returns a non-const pointer to the array of values.
73       * This function is aimed at interoperability with other libraries.
74       * \sa innerIndexPtr(), outerIndexPtr() */
75     inline Scalar* valuePtr() { return derived().valuePtr(); }
76 
77     /** \returns a const pointer to the array of inner indices.
78       * This function is aimed at interoperability with other libraries.
79       * \sa valuePtr(), outerIndexPtr() */
80     inline const StorageIndex* innerIndexPtr() const { return derived().innerIndexPtr(); }
81     /** \returns a non-const pointer to the array of inner indices.
82       * This function is aimed at interoperability with other libraries.
83       * \sa valuePtr(), outerIndexPtr() */
84     inline StorageIndex* innerIndexPtr() { return derived().innerIndexPtr(); }
85 
86     /** \returns a const pointer to the array of the starting positions of the inner vectors.
87       * This function is aimed at interoperability with other libraries.
88       * \warning it returns the null pointer 0 for SparseVector
89       * \sa valuePtr(), innerIndexPtr() */
90     inline const StorageIndex* outerIndexPtr() const { return derived().outerIndexPtr(); }
91     /** \returns a non-const pointer to the array of the starting positions of the inner vectors.
92       * This function is aimed at interoperability with other libraries.
93       * \warning it returns the null pointer 0 for SparseVector
94       * \sa valuePtr(), innerIndexPtr() */
95     inline StorageIndex* outerIndexPtr() { return derived().outerIndexPtr(); }
96 
97     /** \returns a const pointer to the array of the number of non zeros of the inner vectors.
98       * This function is aimed at interoperability with other libraries.
99       * \warning it returns the null pointer 0 in compressed mode */
100     inline const StorageIndex* innerNonZeroPtr() const { return derived().innerNonZeroPtr(); }
101     /** \returns a non-const pointer to the array of the number of non zeros of the inner vectors.
102       * This function is aimed at interoperability with other libraries.
103       * \warning it returns the null pointer 0 in compressed mode */
104     inline StorageIndex* innerNonZeroPtr() { return derived().innerNonZeroPtr(); }
105 
106     /** \returns whether \c *this is in compressed form. */
107     inline bool isCompressed() const { return innerNonZeroPtr()==0; }
108 
109     /** \returns a read-only view of the stored coefficients as a 1D array expression.
110       *
111       * \warning this method is for \b compressed \b storage \b only, and it will trigger an assertion otherwise.
112       *
113       * \sa valuePtr(), isCompressed() */
114     const Map<const Array<Scalar,Dynamic,1> > coeffs() const { eigen_assert(isCompressed()); return Array<Scalar,Dynamic,1>::Map(valuePtr(),nonZeros()); }
115 
116     /** \returns a read-write view of the stored coefficients as a 1D array expression
117       *
118       * \warning this method is for \b compressed \b storage \b only, and it will trigger an assertion otherwise.
119       *
120       * Here is an example:
121       * \include SparseMatrix_coeffs.cpp
122       * and the output is:
123       * \include SparseMatrix_coeffs.out
124       *
125       * \sa valuePtr(), isCompressed() */
126     Map<Array<Scalar,Dynamic,1> > coeffs() { eigen_assert(isCompressed()); return Array<Scalar,Dynamic,1>::Map(valuePtr(),nonZeros()); }
127 
128   protected:
129     /** Default constructor. Do nothing. */
130     SparseCompressedBase() {}
131   private:
132     template<typename OtherDerived> explicit SparseCompressedBase(const SparseCompressedBase<OtherDerived>&);
133 };
134 
135 template<typename Derived>
136 class SparseCompressedBase<Derived>::InnerIterator
137 {
138   public:
139     InnerIterator()
140       : m_values(0), m_indices(0), m_outer(0), m_id(0), m_end(0)
141     {}
142 
143     InnerIterator(const InnerIterator& other)
144       : m_values(other.m_values), m_indices(other.m_indices), m_outer(other.m_outer), m_id(other.m_id), m_end(other.m_end)
145     {}
146 
147     InnerIterator& operator=(const InnerIterator& other)
148     {
149       m_values = other.m_values;
150       m_indices = other.m_indices;
151       const_cast<OuterType&>(m_outer).setValue(other.m_outer.value());
152       m_id = other.m_id;
153       m_end = other.m_end;
154       return *this;
155     }
156 
157     InnerIterator(const SparseCompressedBase& mat, Index outer)
158       : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(outer)
159     {
160       if(Derived::IsVectorAtCompileTime && mat.outerIndexPtr()==0)
161       {
162         m_id = 0;
163         m_end = mat.nonZeros();
164       }
165       else
166       {
167         m_id = mat.outerIndexPtr()[outer];
168         if(mat.isCompressed())
169           m_end = mat.outerIndexPtr()[outer+1];
170         else
171           m_end = m_id + mat.innerNonZeroPtr()[outer];
172       }
173     }
174 
175     explicit InnerIterator(const SparseCompressedBase& mat)
176       : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(0), m_id(0), m_end(mat.nonZeros())
177     {
178       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived);
179     }
180 
181     explicit InnerIterator(const internal::CompressedStorage<Scalar,StorageIndex>& data)
182       : m_values(data.valuePtr()), m_indices(data.indexPtr()), m_outer(0), m_id(0), m_end(data.size())
183     {
184       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived);
185     }
186 
187     inline InnerIterator& operator++() { m_id++; return *this; }
188 
189     inline const Scalar& value() const { return m_values[m_id]; }
190     inline Scalar& valueRef() { return const_cast<Scalar&>(m_values[m_id]); }
191 
192     inline StorageIndex index() const { return m_indices[m_id]; }
193     inline Index outer() const { return m_outer.value(); }
194     inline Index row() const { return IsRowMajor ? m_outer.value() : index(); }
195     inline Index col() const { return IsRowMajor ? index() : m_outer.value(); }
196 
197     inline operator bool() const { return (m_id < m_end); }
198 
199   protected:
200     const Scalar* m_values;
201     const StorageIndex* m_indices;
202     typedef internal::variable_if_dynamic<Index,Derived::IsVectorAtCompileTime?0:Dynamic> OuterType;
203     const OuterType m_outer;
204     Index m_id;
205     Index m_end;
206   private:
207     // If you get here, then you're not using the right InnerIterator type, e.g.:
208     //   SparseMatrix<double,RowMajor> A;
209     //   SparseMatrix<double>::InnerIterator it(A,0);
210     template<typename T> InnerIterator(const SparseMatrixBase<T>&, Index outer);
211 };
212 
213 template<typename Derived>
214 class SparseCompressedBase<Derived>::ReverseInnerIterator
215 {
216   public:
217     ReverseInnerIterator(const SparseCompressedBase& mat, Index outer)
218       : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(outer)
219     {
220       if(Derived::IsVectorAtCompileTime && mat.outerIndexPtr()==0)
221       {
222         m_start = 0;
223         m_id = mat.nonZeros();
224       }
225       else
226       {
227         m_start = mat.outerIndexPtr()[outer];
228         if(mat.isCompressed())
229           m_id = mat.outerIndexPtr()[outer+1];
230         else
231           m_id = m_start + mat.innerNonZeroPtr()[outer];
232       }
233     }
234 
235     explicit ReverseInnerIterator(const SparseCompressedBase& mat)
236       : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(0), m_start(0), m_id(mat.nonZeros())
237     {
238       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived);
239     }
240 
241     explicit ReverseInnerIterator(const internal::CompressedStorage<Scalar,StorageIndex>& data)
242       : m_values(data.valuePtr()), m_indices(data.indexPtr()), m_outer(0), m_start(0), m_id(data.size())
243     {
244       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived);
245     }
246 
247     inline ReverseInnerIterator& operator--() { --m_id; return *this; }
248 
249     inline const Scalar& value() const { return m_values[m_id-1]; }
250     inline Scalar& valueRef() { return const_cast<Scalar&>(m_values[m_id-1]); }
251 
252     inline StorageIndex index() const { return m_indices[m_id-1]; }
253     inline Index outer() const { return m_outer.value(); }
254     inline Index row() const { return IsRowMajor ? m_outer.value() : index(); }
255     inline Index col() const { return IsRowMajor ? index() : m_outer.value(); }
256 
257     inline operator bool() const { return (m_id > m_start); }
258 
259   protected:
260     const Scalar* m_values;
261     const StorageIndex* m_indices;
262     typedef internal::variable_if_dynamic<Index,Derived::IsVectorAtCompileTime?0:Dynamic> OuterType;
263     const OuterType m_outer;
264     Index m_start;
265     Index m_id;
266 };
267 
268 namespace internal {
269 
270 template<typename Derived>
271 struct evaluator<SparseCompressedBase<Derived> >
272   : evaluator_base<Derived>
273 {
274   typedef typename Derived::Scalar Scalar;
275   typedef typename Derived::InnerIterator InnerIterator;
276 
277   enum {
278     CoeffReadCost = NumTraits<Scalar>::ReadCost,
279     Flags = Derived::Flags
280   };
281 
282   evaluator() : m_matrix(0), m_zero(0)
283   {
284     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
285   }
286   explicit evaluator(const Derived &mat) : m_matrix(&mat), m_zero(0)
287   {
288     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
289   }
290 
291   inline Index nonZerosEstimate() const {
292     return m_matrix->nonZeros();
293   }
294 
295   operator Derived&() { return m_matrix->const_cast_derived(); }
296   operator const Derived&() const { return *m_matrix; }
297 
298   typedef typename DenseCoeffsBase<Derived,ReadOnlyAccessors>::CoeffReturnType CoeffReturnType;
299   const Scalar& coeff(Index row, Index col) const
300   {
301     Index p = find(row,col);
302 
303     if(p==Dynamic)
304       return m_zero;
305     else
306       return m_matrix->const_cast_derived().valuePtr()[p];
307   }
308 
309   Scalar& coeffRef(Index row, Index col)
310   {
311     Index p = find(row,col);
312     eigen_assert(p!=Dynamic && "written coefficient does not exist");
313     return m_matrix->const_cast_derived().valuePtr()[p];
314   }
315 
316 protected:
317 
318   Index find(Index row, Index col) const
319   {
320     eigen_internal_assert(row>=0 && row<m_matrix->rows() && col>=0 && col<m_matrix->cols());
321 
322     const Index outer = Derived::IsRowMajor ? row : col;
323     const Index inner = Derived::IsRowMajor ? col : row;
324 
325     Index start = m_matrix->outerIndexPtr()[outer];
326     Index end = m_matrix->isCompressed() ? m_matrix->outerIndexPtr()[outer+1] : m_matrix->outerIndexPtr()[outer] + m_matrix->innerNonZeroPtr()[outer];
327     eigen_assert(end>=start && "you are using a non finalized sparse matrix or written coefficient does not exist");
328     const Index p = std::lower_bound(m_matrix->innerIndexPtr()+start, m_matrix->innerIndexPtr()+end,inner) - m_matrix->innerIndexPtr();
329 
330     return ((p<end) && (m_matrix->innerIndexPtr()[p]==inner)) ? p : Dynamic;
331   }
332 
333   const Derived *m_matrix;
334   const Scalar m_zero;
335 };
336 
337 }
338 
339 } // end namespace Eigen
340 
341 #endif // EIGEN_SPARSE_COMPRESSED_BASE_H
342