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
13namespace 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
35namespace internal {
36template<typename _Scalar, int _Options, typename _StorageIndex>
37struct traits<DynamicSparseMatrix<_Scalar, _Options, _StorageIndex> >
38{
39  typedef _Scalar Scalar;
40  typedef _StorageIndex StorageIndex;
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
55template<typename _Scalar, int _Options, typename _StorageIndex>
56 class  DynamicSparseMatrix
57  : public SparseMatrixBase<DynamicSparseMatrix<_Scalar, _Options, _StorageIndex> >
58{
59    typedef SparseMatrixBase<DynamicSparseMatrix> Base;
60    using Base::convert_index;
61  public:
62    EIGEN_SPARSE_PUBLIC_INTERFACE(DynamicSparseMatrix)
63    // FIXME: why are these operator already alvailable ???
64    // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, +=)
65    // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=)
66    typedef MappedSparseMatrix<Scalar,Flags> Map;
67    using Base::IsRowMajor;
68    using Base::operator=;
69    enum {
70      Options = _Options
71    };
72
73  protected:
74
75    typedef DynamicSparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0), StorageIndex> TransposedSparseMatrix;
76
77    Index m_innerSize;
78    std::vector<internal::CompressedStorage<Scalar,StorageIndex> > m_data;
79
80  public:
81
82    inline Index rows() const { return IsRowMajor ? outerSize() : m_innerSize; }
83    inline Index cols() const { return IsRowMajor ? m_innerSize : outerSize(); }
84    inline Index innerSize() const { return m_innerSize; }
85    inline Index outerSize() const { return convert_index(m_data.size()); }
86    inline Index innerNonZeros(Index j) const { return m_data[j].size(); }
87
88    std::vector<internal::CompressedStorage<Scalar,StorageIndex> >& _data() { return m_data; }
89    const std::vector<internal::CompressedStorage<Scalar,StorageIndex> >& _data() const { return m_data; }
90
91    /** \returns the coefficient value at given position \a row, \a col
92      * This operation involes a log(rho*outer_size) binary search.
93      */
94    inline Scalar coeff(Index row, Index col) const
95    {
96      const Index outer = IsRowMajor ? row : col;
97      const Index inner = IsRowMajor ? col : row;
98      return m_data[outer].at(inner);
99    }
100
101    /** \returns a reference to the coefficient value at given position \a row, \a col
102      * This operation involes a log(rho*outer_size) binary search. If the coefficient does not
103      * exist yet, then a sorted insertion into a sequential buffer is performed.
104      */
105    inline Scalar& coeffRef(Index row, Index col)
106    {
107      const Index outer = IsRowMajor ? row : col;
108      const Index inner = IsRowMajor ? col : row;
109      return m_data[outer].atWithInsertion(inner);
110    }
111
112    class InnerIterator;
113    class ReverseInnerIterator;
114
115    void setZero()
116    {
117      for (Index j=0; j<outerSize(); ++j)
118        m_data[j].clear();
119    }
120
121    /** \returns the number of non zero coefficients */
122    Index nonZeros() const
123    {
124      Index res = 0;
125      for (Index j=0; j<outerSize(); ++j)
126        res += m_data[j].size();
127      return res;
128    }
129
130
131
132    void reserve(Index reserveSize = 1000)
133    {
134      if (outerSize()>0)
135      {
136        Index reserveSizePerVector = (std::max)(reserveSize/outerSize(),Index(4));
137        for (Index j=0; j<outerSize(); ++j)
138        {
139          m_data[j].reserve(reserveSizePerVector);
140        }
141      }
142    }
143
144    /** Does nothing: provided for compatibility with SparseMatrix */
145    inline void startVec(Index /*outer*/) {}
146
147    /** \returns a reference to the non zero coefficient at position \a row, \a col assuming that:
148      * - the nonzero does not already exist
149      * - the new coefficient is the last one of the given inner vector.
150      *
151      * \sa insert, insertBackByOuterInner */
152    inline Scalar& insertBack(Index row, Index col)
153    {
154      return insertBackByOuterInner(IsRowMajor?row:col, IsRowMajor?col:row);
155    }
156
157    /** \sa insertBack */
158    inline Scalar& insertBackByOuterInner(Index outer, Index inner)
159    {
160      eigen_assert(outer<Index(m_data.size()) && inner<m_innerSize && "out of range");
161      eigen_assert(((m_data[outer].size()==0) || (m_data[outer].index(m_data[outer].size()-1)<inner))
162                && "wrong sorted insertion");
163      m_data[outer].append(0, inner);
164      return m_data[outer].value(m_data[outer].size()-1);
165    }
166
167    inline Scalar& insert(Index row, Index col)
168    {
169      const Index outer = IsRowMajor ? row : col;
170      const Index inner = IsRowMajor ? col : row;
171
172      Index startId = 0;
173      Index id = static_cast<Index>(m_data[outer].size()) - 1;
174      m_data[outer].resize(id+2,1);
175
176      while ( (id >= startId) && (m_data[outer].index(id) > inner) )
177      {
178        m_data[outer].index(id+1) = m_data[outer].index(id);
179        m_data[outer].value(id+1) = m_data[outer].value(id);
180        --id;
181      }
182      m_data[outer].index(id+1) = inner;
183      m_data[outer].value(id+1) = 0;
184      return m_data[outer].value(id+1);
185    }
186
187    /** Does nothing: provided for compatibility with SparseMatrix */
188    inline void finalize() {}
189
190    /** Suppress all nonzeros which are smaller than \a reference under the tolerence \a epsilon */
191    void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
192    {
193      for (Index j=0; j<outerSize(); ++j)
194        m_data[j].prune(reference,epsilon);
195    }
196
197    /** Resize the matrix without preserving the data (the matrix is set to zero)
198      */
199    void resize(Index rows, Index cols)
200    {
201      const Index outerSize = IsRowMajor ? rows : cols;
202      m_innerSize = convert_index(IsRowMajor ? cols : rows);
203      setZero();
204      if (Index(m_data.size()) != outerSize)
205      {
206        m_data.resize(outerSize);
207      }
208    }
209
210    void resizeAndKeepData(Index rows, Index cols)
211    {
212      const Index outerSize = IsRowMajor ? rows : cols;
213      const Index innerSize = IsRowMajor ? cols : rows;
214      if (m_innerSize>innerSize)
215      {
216        // remove all coefficients with innerCoord>=innerSize
217        // TODO
218        //std::cerr << "not implemented yet\n";
219        exit(2);
220      }
221      if (m_data.size() != outerSize)
222      {
223        m_data.resize(outerSize);
224      }
225    }
226
227    /** The class DynamicSparseMatrix is deprectaed */
228    EIGEN_DEPRECATED inline DynamicSparseMatrix()
229      : m_innerSize(0), m_data(0)
230    {
231      eigen_assert(innerSize()==0 && outerSize()==0);
232    }
233
234    /** The class DynamicSparseMatrix is deprectaed */
235    EIGEN_DEPRECATED inline DynamicSparseMatrix(Index rows, Index cols)
236      : m_innerSize(0)
237    {
238      resize(rows, cols);
239    }
240
241    /** The class DynamicSparseMatrix is deprectaed */
242    template<typename OtherDerived>
243    EIGEN_DEPRECATED explicit inline DynamicSparseMatrix(const SparseMatrixBase<OtherDerived>& other)
244      : m_innerSize(0)
245    {
246    Base::operator=(other.derived());
247    }
248
249    inline DynamicSparseMatrix(const DynamicSparseMatrix& other)
250      : Base(), m_innerSize(0)
251    {
252      *this = other.derived();
253    }
254
255    inline void swap(DynamicSparseMatrix& other)
256    {
257      //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
258      std::swap(m_innerSize, other.m_innerSize);
259      //std::swap(m_outerSize, other.m_outerSize);
260      m_data.swap(other.m_data);
261    }
262
263    inline DynamicSparseMatrix& operator=(const DynamicSparseMatrix& other)
264    {
265      if (other.isRValue())
266      {
267        swap(other.const_cast_derived());
268      }
269      else
270      {
271        resize(other.rows(), other.cols());
272        m_data = other.m_data;
273      }
274      return *this;
275    }
276
277    /** Destructor */
278    inline ~DynamicSparseMatrix() {}
279
280  public:
281
282    /** \deprecated
283      * Set the matrix to zero and reserve the memory for \a reserveSize nonzero coefficients. */
284    EIGEN_DEPRECATED void startFill(Index reserveSize = 1000)
285    {
286      setZero();
287      reserve(reserveSize);
288    }
289
290    /** \deprecated use insert()
291      * inserts a nonzero coefficient at given coordinates \a row, \a col and returns its reference assuming that:
292      *  1 - the coefficient does not exist yet
293      *  2 - this the coefficient with greater inner coordinate for the given outer coordinate.
294      * In other words, assuming \c *this is column-major, then there must not exists any nonzero coefficient of coordinates
295      * \c i \c x \a col such that \c i >= \a row. Otherwise the matrix is invalid.
296      *
297      * \see fillrand(), coeffRef()
298      */
299    EIGEN_DEPRECATED Scalar& fill(Index row, Index col)
300    {
301      const Index outer = IsRowMajor ? row : col;
302      const Index inner = IsRowMajor ? col : row;
303      return insertBack(outer,inner);
304    }
305
306    /** \deprecated use insert()
307      * Like fill() but with random inner coordinates.
308      * Compared to the generic coeffRef(), the unique limitation is that we assume
309      * the coefficient does not exist yet.
310      */
311    EIGEN_DEPRECATED Scalar& fillrand(Index row, Index col)
312    {
313      return insert(row,col);
314    }
315
316    /** \deprecated use finalize()
317      * Does nothing. Provided for compatibility with SparseMatrix. */
318    EIGEN_DEPRECATED void endFill() {}
319
320#   ifdef EIGEN_DYNAMICSPARSEMATRIX_PLUGIN
321#     include EIGEN_DYNAMICSPARSEMATRIX_PLUGIN
322#   endif
323 };
324
325template<typename Scalar, int _Options, typename _StorageIndex>
326class DynamicSparseMatrix<Scalar,_Options,_StorageIndex>::InnerIterator : public SparseVector<Scalar,_Options,_StorageIndex>::InnerIterator
327{
328    typedef typename SparseVector<Scalar,_Options,_StorageIndex>::InnerIterator Base;
329  public:
330    InnerIterator(const DynamicSparseMatrix& mat, Index outer)
331      : Base(mat.m_data[outer]), m_outer(outer)
332    {}
333
334    inline Index row() const { return IsRowMajor ? m_outer : Base::index(); }
335    inline Index col() const { return IsRowMajor ? Base::index() : m_outer; }
336    inline Index outer() const { return m_outer; }
337
338  protected:
339    const Index m_outer;
340};
341
342template<typename Scalar, int _Options, typename _StorageIndex>
343class DynamicSparseMatrix<Scalar,_Options,_StorageIndex>::ReverseInnerIterator : public SparseVector<Scalar,_Options,_StorageIndex>::ReverseInnerIterator
344{
345    typedef typename SparseVector<Scalar,_Options,_StorageIndex>::ReverseInnerIterator Base;
346  public:
347    ReverseInnerIterator(const DynamicSparseMatrix& mat, Index outer)
348      : Base(mat.m_data[outer]), m_outer(outer)
349    {}
350
351    inline Index row() const { return IsRowMajor ? m_outer : Base::index(); }
352    inline Index col() const { return IsRowMajor ? Base::index() : m_outer; }
353    inline Index outer() const { return m_outer; }
354
355  protected:
356    const Index m_outer;
357};
358
359namespace internal {
360
361template<typename _Scalar, int _Options, typename _StorageIndex>
362struct evaluator<DynamicSparseMatrix<_Scalar,_Options,_StorageIndex> >
363  : evaluator_base<DynamicSparseMatrix<_Scalar,_Options,_StorageIndex> >
364{
365  typedef _Scalar Scalar;
366  typedef DynamicSparseMatrix<_Scalar,_Options,_StorageIndex> SparseMatrixType;
367  typedef typename SparseMatrixType::InnerIterator InnerIterator;
368  typedef typename SparseMatrixType::ReverseInnerIterator ReverseInnerIterator;
369
370  enum {
371    CoeffReadCost = NumTraits<_Scalar>::ReadCost,
372    Flags = SparseMatrixType::Flags
373  };
374
375  evaluator() : m_matrix(0) {}
376  evaluator(const SparseMatrixType &mat) : m_matrix(&mat) {}
377
378  operator SparseMatrixType&() { return m_matrix->const_cast_derived(); }
379  operator const SparseMatrixType&() const { return *m_matrix; }
380
381  Scalar coeff(Index row, Index col) const { return m_matrix->coeff(row,col); }
382
383  Index nonZerosEstimate() const { return m_matrix->nonZeros(); }
384
385  const SparseMatrixType *m_matrix;
386};
387
388}
389
390} // end namespace Eigen
391
392#endif // EIGEN_DYNAMIC_SPARSEMATRIX_H
393