Home | History | Annotate | Download | only in SparseExtra
      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 (at) 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 _StorageIndex>
     37 struct 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 
     55 template<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 
    325 template<typename Scalar, int _Options, typename _StorageIndex>
    326 class 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 
    342 template<typename Scalar, int _Options, typename _StorageIndex>
    343 class 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 
    359 namespace internal {
    360 
    361 template<typename _Scalar, int _Options, typename _StorageIndex>
    362 struct 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