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_DYNAMIC_SPARSEMATRIX_H
11 #define EIGEN_DYNAMIC_SPARSEMATRIX_H
12 
13 namespace Eigen {
14 
15 /** \deprecated use a SparseMatrix in an uncompressed mode
16   *
17   * \class DynamicSparseMatrix
18   *
19   * \brief A sparse matrix class designed for matrix assembly purpose
20   *
21   * \param _Scalar the scalar type, i.e. the type of the coefficients
22   *
23   * Unlike SparseMatrix, this class provides a much higher degree of flexibility. In particular, it allows
24   * random read/write accesses in log(rho*outer_size) where \c rho is the probability that a coefficient is
25   * nonzero and outer_size is the number of columns if the matrix is column-major and the number of rows
26   * otherwise.
27   *
28   * Internally, the data are stored as a std::vector of compressed vector. The performances of random writes might
29   * decrease as the number of nonzeros per inner-vector increase. In practice, we observed very good performance
30   * till about 100 nonzeros/vector, and the performance remains relatively good till 500 nonzeros/vectors.
31   *
32   * \see SparseMatrix
33   */
34 
35 namespace internal {
36 template<typename _Scalar, int _Options, typename _Index>
37 struct traits<DynamicSparseMatrix<_Scalar, _Options, _Index> >
38 {
39   typedef _Scalar Scalar;
40   typedef _Index Index;
41   typedef Sparse StorageKind;
42   typedef MatrixXpr XprKind;
43   enum {
44     RowsAtCompileTime = Dynamic,
45     ColsAtCompileTime = Dynamic,
46     MaxRowsAtCompileTime = Dynamic,
47     MaxColsAtCompileTime = Dynamic,
48     Flags = _Options | NestByRefBit | LvalueBit,
49     CoeffReadCost = NumTraits<Scalar>::ReadCost,
50     SupportedAccessPatterns = OuterRandomAccessPattern
51   };
52 };
53 }
54 
55 template<typename _Scalar, int _Options, typename _Index>
56  class  DynamicSparseMatrix
57   : public SparseMatrixBase<DynamicSparseMatrix<_Scalar, _Options, _Index> >
58 {
59   public:
60     EIGEN_SPARSE_PUBLIC_INTERFACE(DynamicSparseMatrix)
61     // FIXME: why are these operator already alvailable ???
62     // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, +=)
63     // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=)
64     typedef MappedSparseMatrix<Scalar,Flags> Map;
65     using Base::IsRowMajor;
66     using Base::operator=;
67     enum {
68       Options = _Options
69     };
70 
71   protected:
72 
73     typedef DynamicSparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix;
74 
75     Index m_innerSize;
76     std::vector<internal::CompressedStorage<Scalar,Index> > m_data;
77 
78   public:
79 
80     inline Index rows() const { return IsRowMajor ? outerSize() : m_innerSize; }
81     inline Index cols() const { return IsRowMajor ? m_innerSize : outerSize(); }
82     inline Index innerSize() const { return m_innerSize; }
83     inline Index outerSize() const { return static_cast<Index>(m_data.size()); }
84     inline Index innerNonZeros(Index j) const { return m_data[j].size(); }
85 
86     std::vector<internal::CompressedStorage<Scalar,Index> >& _data() { return m_data; }
87     const std::vector<internal::CompressedStorage<Scalar,Index> >& _data() const { return m_data; }
88 
89     /** \returns the coefficient value at given position \a row, \a col
90       * This operation involes a log(rho*outer_size) binary search.
91       */
92     inline Scalar coeff(Index row, Index col) const
93     {
94       const Index outer = IsRowMajor ? row : col;
95       const Index inner = IsRowMajor ? col : row;
96       return m_data[outer].at(inner);
97     }
98 
99     /** \returns a reference to the coefficient value at given position \a row, \a col
100       * This operation involes a log(rho*outer_size) binary search. If the coefficient does not
101       * exist yet, then a sorted insertion into a sequential buffer is performed.
102       */
103     inline Scalar& coeffRef(Index row, Index col)
104     {
105       const Index outer = IsRowMajor ? row : col;
106       const Index inner = IsRowMajor ? col : row;
107       return m_data[outer].atWithInsertion(inner);
108     }
109 
110     class InnerIterator;
111     class ReverseInnerIterator;
112 
113     void setZero()
114     {
115       for (Index j=0; j<outerSize(); ++j)
116         m_data[j].clear();
117     }
118 
119     /** \returns the number of non zero coefficients */
120     Index nonZeros() const
121     {
122       Index res = 0;
123       for (Index j=0; j<outerSize(); ++j)
124         res += static_cast<Index>(m_data[j].size());
125       return res;
126     }
127 
128 
129 
130     void reserve(Index reserveSize = 1000)
131     {
132       if (outerSize()>0)
133       {
134         Index reserveSizePerVector = (std::max)(reserveSize/outerSize(),Index(4));
135         for (Index j=0; j<outerSize(); ++j)
136         {
137           m_data[j].reserve(reserveSizePerVector);
138         }
139       }
140     }
141 
142     /** Does nothing: provided for compatibility with SparseMatrix */
143     inline void startVec(Index /*outer*/) {}
144 
145     /** \returns a reference to the non zero coefficient at position \a row, \a col assuming that:
146       * - the nonzero does not already exist
147       * - the new coefficient is the last one of the given inner vector.
148       *
149       * \sa insert, insertBackByOuterInner */
150     inline Scalar& insertBack(Index row, Index col)
151     {
152       return insertBackByOuterInner(IsRowMajor?row:col, IsRowMajor?col:row);
153     }
154 
155     /** \sa insertBack */
156     inline Scalar& insertBackByOuterInner(Index outer, Index inner)
157     {
158       eigen_assert(outer<Index(m_data.size()) && inner<m_innerSize && "out of range");
159       eigen_assert(((m_data[outer].size()==0) || (m_data[outer].index(m_data[outer].size()-1)<inner))
160                 && "wrong sorted insertion");
161       m_data[outer].append(0, inner);
162       return m_data[outer].value(m_data[outer].size()-1);
163     }
164 
165     inline Scalar& insert(Index row, Index col)
166     {
167       const Index outer = IsRowMajor ? row : col;
168       const Index inner = IsRowMajor ? col : row;
169 
170       Index startId = 0;
171       Index id = static_cast<Index>(m_data[outer].size()) - 1;
172       m_data[outer].resize(id+2,1);
173 
174       while ( (id >= startId) && (m_data[outer].index(id) > inner) )
175       {
176         m_data[outer].index(id+1) = m_data[outer].index(id);
177         m_data[outer].value(id+1) = m_data[outer].value(id);
178         --id;
179       }
180       m_data[outer].index(id+1) = inner;
181       m_data[outer].value(id+1) = 0;
182       return m_data[outer].value(id+1);
183     }
184 
185     /** Does nothing: provided for compatibility with SparseMatrix */
186     inline void finalize() {}
187 
188     /** Suppress all nonzeros which are smaller than \a reference under the tolerence \a epsilon */
189     void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
190     {
191       for (Index j=0; j<outerSize(); ++j)
192         m_data[j].prune(reference,epsilon);
193     }
194 
195     /** Resize the matrix without preserving the data (the matrix is set to zero)
196       */
197     void resize(Index rows, Index cols)
198     {
199       const Index outerSize = IsRowMajor ? rows : cols;
200       m_innerSize = IsRowMajor ? cols : rows;
201       setZero();
202       if (Index(m_data.size()) != outerSize)
203       {
204         m_data.resize(outerSize);
205       }
206     }
207 
208     void resizeAndKeepData(Index rows, Index cols)
209     {
210       const Index outerSize = IsRowMajor ? rows : cols;
211       const Index innerSize = IsRowMajor ? cols : rows;
212       if (m_innerSize>innerSize)
213       {
214         // remove all coefficients with innerCoord>=innerSize
215         // TODO
216         //std::cerr << "not implemented yet\n";
217         exit(2);
218       }
219       if (m_data.size() != outerSize)
220       {
221         m_data.resize(outerSize);
222       }
223     }
224 
225     /** The class DynamicSparseMatrix is deprectaed */
226     EIGEN_DEPRECATED inline DynamicSparseMatrix()
227       : m_innerSize(0), m_data(0)
228     {
229       eigen_assert(innerSize()==0 && outerSize()==0);
230     }
231 
232     /** The class DynamicSparseMatrix is deprectaed */
233     EIGEN_DEPRECATED inline DynamicSparseMatrix(Index rows, Index cols)
234       : m_innerSize(0)
235     {
236       resize(rows, cols);
237     }
238 
239     /** The class DynamicSparseMatrix is deprectaed */
240     template<typename OtherDerived>
241     EIGEN_DEPRECATED explicit inline DynamicSparseMatrix(const SparseMatrixBase<OtherDerived>& other)
242       : m_innerSize(0)
243     {
244     Base::operator=(other.derived());
245     }
246 
247     inline DynamicSparseMatrix(const DynamicSparseMatrix& other)
248       : Base(), m_innerSize(0)
249     {
250       *this = other.derived();
251     }
252 
253     inline void swap(DynamicSparseMatrix& other)
254     {
255       //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
256       std::swap(m_innerSize, other.m_innerSize);
257       //std::swap(m_outerSize, other.m_outerSize);
258       m_data.swap(other.m_data);
259     }
260 
261     inline DynamicSparseMatrix& operator=(const DynamicSparseMatrix& other)
262     {
263       if (other.isRValue())
264       {
265         swap(other.const_cast_derived());
266       }
267       else
268       {
269         resize(other.rows(), other.cols());
270         m_data = other.m_data;
271       }
272       return *this;
273     }
274 
275     /** Destructor */
276     inline ~DynamicSparseMatrix() {}
277 
278   public:
279 
280     /** \deprecated
281       * Set the matrix to zero and reserve the memory for \a reserveSize nonzero coefficients. */
282     EIGEN_DEPRECATED void startFill(Index reserveSize = 1000)
283     {
284       setZero();
285       reserve(reserveSize);
286     }
287 
288     /** \deprecated use insert()
289       * inserts a nonzero coefficient at given coordinates \a row, \a col and returns its reference assuming that:
290       *  1 - the coefficient does not exist yet
291       *  2 - this the coefficient with greater inner coordinate for the given outer coordinate.
292       * In other words, assuming \c *this is column-major, then there must not exists any nonzero coefficient of coordinates
293       * \c i \c x \a col such that \c i >= \a row. Otherwise the matrix is invalid.
294       *
295       * \see fillrand(), coeffRef()
296       */
297     EIGEN_DEPRECATED Scalar& fill(Index row, Index col)
298     {
299       const Index outer = IsRowMajor ? row : col;
300       const Index inner = IsRowMajor ? col : row;
301       return insertBack(outer,inner);
302     }
303 
304     /** \deprecated use insert()
305       * Like fill() but with random inner coordinates.
306       * Compared to the generic coeffRef(), the unique limitation is that we assume
307       * the coefficient does not exist yet.
308       */
309     EIGEN_DEPRECATED Scalar& fillrand(Index row, Index col)
310     {
311       return insert(row,col);
312     }
313 
314     /** \deprecated use finalize()
315       * Does nothing. Provided for compatibility with SparseMatrix. */
316     EIGEN_DEPRECATED void endFill() {}
317 
318 #   ifdef EIGEN_DYNAMICSPARSEMATRIX_PLUGIN
319 #     include EIGEN_DYNAMICSPARSEMATRIX_PLUGIN
320 #   endif
321  };
322 
323 template<typename Scalar, int _Options, typename _Index>
324 class DynamicSparseMatrix<Scalar,_Options,_Index>::InnerIterator : public SparseVector<Scalar,_Options,_Index>::InnerIterator
325 {
326     typedef typename SparseVector<Scalar,_Options,_Index>::InnerIterator Base;
327   public:
328     InnerIterator(const DynamicSparseMatrix& mat, Index outer)
329       : Base(mat.m_data[outer]), m_outer(outer)
330     {}
331 
332     inline Index row() const { return IsRowMajor ? m_outer : Base::index(); }
333     inline Index col() const { return IsRowMajor ? Base::index() : m_outer; }
334 
335   protected:
336     const Index m_outer;
337 };
338 
339 template<typename Scalar, int _Options, typename _Index>
340 class DynamicSparseMatrix<Scalar,_Options,_Index>::ReverseInnerIterator : public SparseVector<Scalar,_Options,_Index>::ReverseInnerIterator
341 {
342     typedef typename SparseVector<Scalar,_Options,_Index>::ReverseInnerIterator Base;
343   public:
344     ReverseInnerIterator(const DynamicSparseMatrix& mat, Index outer)
345       : Base(mat.m_data[outer]), m_outer(outer)
346     {}
347 
348     inline Index row() const { return IsRowMajor ? m_outer : Base::index(); }
349     inline Index col() const { return IsRowMajor ? Base::index() : m_outer; }
350 
351   protected:
352     const Index m_outer;
353 };
354 
355 } // end namespace Eigen
356 
357 #endif // EIGEN_DYNAMIC_SPARSEMATRIX_H
358