Home | History | Annotate | Download | only in SparseCore 
      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_SPARSE_BLOCK_H
     11 #define EIGEN_SPARSE_BLOCK_H
     12 
     13 namespace Eigen {
     14 
     15 template<typename XprType, int BlockRows, int BlockCols>
     16 class BlockImpl<XprType,BlockRows,BlockCols,true,Sparse>
     17   : public SparseMatrixBase<Block<XprType,BlockRows,BlockCols,true> >
     18 {
     19     typedef typename internal::remove_all<typename XprType::Nested>::type _MatrixTypeNested;
     20     typedef Block<XprType, BlockRows, BlockCols, true> BlockType;
     21 public:
     22     enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
     23 protected:
     24     enum { OuterSize = IsRowMajor ? BlockRows : BlockCols };
     25 public:
     26     EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType)
     27 
     28     class InnerIterator: public XprType::InnerIterator
     29     {
     30         typedef typename BlockImpl::Index Index;
     31       public:
     32         inline InnerIterator(const BlockType& xpr, Index outer)
     33           : XprType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
     34         {}
     35         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
     36         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
     37       protected:
     38         Index m_outer;
     39     };
     40     class ReverseInnerIterator: public XprType::ReverseInnerIterator
     41     {
     42         typedef typename BlockImpl::Index Index;
     43       public:
     44         inline ReverseInnerIterator(const BlockType& xpr, Index outer)
     45           : XprType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
     46         {}
     47         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
     48         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
     49       protected:
     50         Index m_outer;
     51     };
     52 
     53     inline BlockImpl(const XprType& xpr, int i)
     54       : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize)
     55     {}
     56 
     57     inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols)
     58       : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols)
     59     {}
     60 
     61     EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
     62     EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
     63 
     64   protected:
     65 
     66     typename XprType::Nested m_matrix;
     67     Index m_outerStart;
     68     const internal::variable_if_dynamic<Index, OuterSize> m_outerSize;
     69 
     70     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl)
     71 };
     72 
     73 
     74 /***************************************************************************
     75 * specialisation for SparseMatrix
     76 ***************************************************************************/
     77 
     78 template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols>
     79 class BlockImpl<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true,Sparse>
     80   : public SparseMatrixBase<Block<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true> >
     81 {
     82     typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType;
     83     typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested;
     84     typedef Block<SparseMatrixType, BlockRows, BlockCols, true> BlockType;
     85 public:
     86     enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
     87     EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType)
     88 protected:
     89     enum { OuterSize = IsRowMajor ? BlockRows : BlockCols };
     90 public:
     91 
     92     class InnerIterator: public SparseMatrixType::InnerIterator
     93     {
     94       public:
     95         inline InnerIterator(const BlockType& xpr, Index outer)
     96           : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
     97         {}
     98         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
     99         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
    100       protected:
    101         Index m_outer;
    102     };
    103     class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator
    104     {
    105       public:
    106         inline ReverseInnerIterator(const BlockType& xpr, Index outer)
    107           : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
    108         {}
    109         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
    110         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
    111       protected:
    112         Index m_outer;
    113     };
    114 
    115     inline BlockImpl(const SparseMatrixType& xpr, int i)
    116       : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize)
    117     {}
    118 
    119     inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols)
    120       : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols)
    121     {}
    122 
    123     template<typename OtherDerived>
    124     inline BlockType& operator=(const SparseMatrixBase<OtherDerived>& other)
    125     {
    126       typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _NestedMatrixType;
    127       _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);;
    128       // This assignement is slow if this vector set is not empty
    129       // and/or it is not at the end of the nonzeros of the underlying matrix.
    130 
    131       // 1 - eval to a temporary to avoid transposition and/or aliasing issues
    132       SparseMatrix<Scalar, IsRowMajor ? RowMajor : ColMajor, Index> tmp(other);
    133 
    134       // 2 - let's check whether there is enough allocated memory
    135       Index nnz           = tmp.nonZeros();
    136       Index start         = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block
    137       Index end           = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]; // ending posiiton of the current block
    138       Index block_size    = end - start;                                                // available room in the current block
    139       Index tail_size     = m_matrix.outerIndexPtr()[m_matrix.outerSize()] - end;
    140 
    141       Index free_size     = m_matrix.isCompressed()
    142                           ? Index(matrix.data().allocatedSize()) + block_size
    143                           : block_size;
    144 
    145       if(nnz>free_size)
    146       {
    147         // realloc manually to reduce copies
    148         typename SparseMatrixType::Storage newdata(m_matrix.data().allocatedSize() - block_size + nnz);
    149 
    150         std::memcpy(&newdata.value(0), &m_matrix.data().value(0), start*sizeof(Scalar));
    151         std::memcpy(&newdata.index(0), &m_matrix.data().index(0), start*sizeof(Index));
    152 
    153         std::memcpy(&newdata.value(start), &tmp.data().value(0), nnz*sizeof(Scalar));
    154         std::memcpy(&newdata.index(start), &tmp.data().index(0), nnz*sizeof(Index));
    155 
    156         std::memcpy(&newdata.value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar));
    157         std::memcpy(&newdata.index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index));
    158 
    159         newdata.resize(m_matrix.outerIndexPtr()[m_matrix.outerSize()] - block_size + nnz);
    160 
    161         matrix.data().swap(newdata);
    162       }
    163       else
    164       {
    165         // no need to realloc, simply copy the tail at its respective position and insert tmp
    166         matrix.data().resize(start + nnz + tail_size);
    167 
    168         std::memmove(&matrix.data().value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar));
    169         std::memmove(&matrix.data().index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index));
    170 
    171         std::memcpy(&matrix.data().value(start), &tmp.data().value(0), nnz*sizeof(Scalar));
    172         std::memcpy(&matrix.data().index(start), &tmp.data().index(0), nnz*sizeof(Index));
    173       }
    174 
    175       // update innerNonZeros
    176       if(!m_matrix.isCompressed())
    177         for(Index j=0; j<m_outerSize.value(); ++j)
    178           matrix.innerNonZeroPtr()[m_outerStart+j] = tmp.innerVector(j).nonZeros();
    179 
    180       // update outer index pointers
    181       Index p = start;
    182       for(Index k=0; k<m_outerSize.value(); ++k)
    183       {
    184         matrix.outerIndexPtr()[m_outerStart+k] = p;
    185         p += tmp.innerVector(k).nonZeros();
    186       }
    187       std::ptrdiff_t offset = nnz - block_size;
    188       for(Index k = m_outerStart + m_outerSize.value(); k<=matrix.outerSize(); ++k)
    189       {
    190         matrix.outerIndexPtr()[k] += offset;
    191       }
    192 
    193       return derived();
    194     }
    195 
    196     inline BlockType& operator=(const BlockType& other)
    197     {
    198       return operator=<BlockType>(other);
    199     }
    200 
    201     inline const Scalar* valuePtr() const
    202     { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
    203     inline Scalar* valuePtr()
    204     { return m_matrix.const_cast_derived().valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
    205 
    206     inline const Index* innerIndexPtr() const
    207     { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
    208     inline Index* innerIndexPtr()
    209     { return m_matrix.const_cast_derived().innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
    210 
    211     inline const Index* outerIndexPtr() const
    212     { return m_matrix.outerIndexPtr() + m_outerStart; }
    213     inline Index* outerIndexPtr()
    214     { return m_matrix.const_cast_derived().outerIndexPtr() + m_outerStart; }
    215 
    216     Index nonZeros() const
    217     {
    218       if(m_matrix.isCompressed())
    219         return  std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()])
    220               - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]);
    221       else if(m_outerSize.value()==0)
    222         return 0;
    223       else
    224         return Map<const Matrix<Index,OuterSize,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum();
    225     }
    226 
    227     const Scalar& lastCoeff() const
    228     {
    229       EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl);
    230       eigen_assert(nonZeros()>0);
    231       if(m_matrix.isCompressed())
    232         return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1];
    233       else
    234         return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1];
    235     }
    236 
    237     EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
    238     EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
    239 
    240   protected:
    241 
    242     typename SparseMatrixType::Nested m_matrix;
    243     Index m_outerStart;
    244     const internal::variable_if_dynamic<Index, OuterSize> m_outerSize;
    245 
    246 };
    247 
    248 //----------
    249 
    250 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
    251   * is col-major (resp. row-major).
    252   */
    253 template<typename Derived>
    254 typename SparseMatrixBase<Derived>::InnerVectorReturnType SparseMatrixBase<Derived>::innerVector(Index outer)
    255 { return InnerVectorReturnType(derived(), outer); }
    256 
    257 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
    258   * is col-major (resp. row-major). Read-only.
    259   */
    260 template<typename Derived>
    261 const typename SparseMatrixBase<Derived>::ConstInnerVectorReturnType SparseMatrixBase<Derived>::innerVector(Index outer) const
    262 { return ConstInnerVectorReturnType(derived(), outer); }
    263 
    264 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
    265   * is col-major (resp. row-major).
    266   */
    267 template<typename Derived>
    268 Block<Derived,Dynamic,Dynamic,true> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize)
    269 {
    270   return Block<Derived,Dynamic,Dynamic,true>(derived(),
    271                                              IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart,
    272                                              IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize);
    273 
    274 }
    275 
    276 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
    277   * is col-major (resp. row-major). Read-only.
    278   */
    279 template<typename Derived>
    280 const Block<const Derived,Dynamic,Dynamic,true> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize) const
    281 {
    282   return Block<const Derived,Dynamic,Dynamic,true>(derived(),
    283                                                   IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart,
    284                                                   IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize);
    285 
    286 }
    287 
    288 /** Generic implementation of sparse Block expression.
    289   * Real-only.
    290   */
    291 template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
    292 class BlockImpl<XprType,BlockRows,BlockCols,InnerPanel,Sparse>
    293   : public SparseMatrixBase<Block<XprType,BlockRows,BlockCols,InnerPanel> >, internal::no_assignment_operator
    294 {
    295   typedef typename internal::remove_all<typename XprType::Nested>::type _MatrixTypeNested;
    296   typedef Block<XprType, BlockRows, BlockCols, InnerPanel> BlockType;
    297 public:
    298     enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
    299     EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType)
    300 
    301     /** Column or Row constructor
    302       */
    303     inline BlockImpl(const XprType& xpr, int i)
    304       : m_matrix(xpr),
    305         m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0),
    306         m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0),
    307         m_blockRows(xpr.rows()),
    308         m_blockCols(xpr.cols())
    309     {}
    310 
    311     /** Dynamic-size constructor
    312       */
    313     inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols)
    314       : m_matrix(xpr), m_startRow(startRow), m_startCol(startCol), m_blockRows(blockRows), m_blockCols(blockCols)
    315     {}
    316 
    317     inline int rows() const { return m_blockRows.value(); }
    318     inline int cols() const { return m_blockCols.value(); }
    319 
    320     inline Scalar& coeffRef(int row, int col)
    321     {
    322       return m_matrix.const_cast_derived()
    323                .coeffRef(row + m_startRow.value(), col + m_startCol.value());
    324     }
    325 
    326     inline const Scalar coeff(int row, int col) const
    327     {
    328       return m_matrix.coeff(row + m_startRow.value(), col + m_startCol.value());
    329     }
    330 
    331     inline Scalar& coeffRef(int index)
    332     {
    333       return m_matrix.const_cast_derived()
    334              .coeffRef(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index),
    335                        m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0));
    336     }
    337 
    338     inline const Scalar coeff(int index) const
    339     {
    340       return m_matrix
    341              .coeff(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index),
    342                     m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0));
    343     }
    344 
    345     inline const _MatrixTypeNested& nestedExpression() const { return m_matrix; }
    346 
    347     class InnerIterator : public _MatrixTypeNested::InnerIterator
    348     {
    349       typedef typename _MatrixTypeNested::InnerIterator Base;
    350       const BlockType& m_block;
    351       Index m_end;
    352     public:
    353 
    354       EIGEN_STRONG_INLINE InnerIterator(const BlockType& block, Index outer)
    355         : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())),
    356           m_block(block),
    357           m_end(IsRowMajor ? block.m_startCol.value()+block.m_blockCols.value() : block.m_startRow.value()+block.m_blockRows.value())
    358       {
    359         while( (Base::operator bool()) && (Base::index() < (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value())) )
    360           Base::operator++();
    361       }
    362 
    363       inline Index index()  const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); }
    364       inline Index outer()  const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); }
    365       inline Index row()    const { return Base::row()   - m_block.m_startRow.value(); }
    366       inline Index col()    const { return Base::col()   - m_block.m_startCol.value(); }
    367 
    368       inline operator bool() const { return Base::operator bool() && Base::index() < m_end; }
    369     };
    370     class ReverseInnerIterator : public _MatrixTypeNested::ReverseInnerIterator
    371     {
    372       typedef typename _MatrixTypeNested::ReverseInnerIterator Base;
    373       const BlockType& m_block;
    374       Index m_begin;
    375     public:
    376 
    377       EIGEN_STRONG_INLINE ReverseInnerIterator(const BlockType& block, Index outer)
    378         : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())),
    379           m_block(block),
    380           m_begin(IsRowMajor ? block.m_startCol.value() : block.m_startRow.value())
    381       {
    382         while( (Base::operator bool()) && (Base::index() >= (IsRowMajor ? m_block.m_startCol.value()+block.m_blockCols.value() : m_block.m_startRow.value()+block.m_blockRows.value())) )
    383           Base::operator--();
    384       }
    385 
    386       inline Index index()  const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); }
    387       inline Index outer()  const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); }
    388       inline Index row()    const { return Base::row()   - m_block.m_startRow.value(); }
    389       inline Index col()    const { return Base::col()   - m_block.m_startCol.value(); }
    390 
    391       inline operator bool() const { return Base::operator bool() && Base::index() >= m_begin; }
    392     };
    393   protected:
    394     friend class InnerIterator;
    395     friend class ReverseInnerIterator;
    396 
    397     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl)
    398 
    399     typename XprType::Nested m_matrix;
    400     const internal::variable_if_dynamic<Index, XprType::RowsAtCompileTime == 1 ? 0 : Dynamic> m_startRow;
    401     const internal::variable_if_dynamic<Index, XprType::ColsAtCompileTime == 1 ? 0 : Dynamic> m_startCol;
    402     const internal::variable_if_dynamic<Index, RowsAtCompileTime> m_blockRows;
    403     const internal::variable_if_dynamic<Index, ColsAtCompileTime> m_blockCols;
    404 
    405 };
    406 
    407 } // end namespace Eigen
    408 
    409 #endif // EIGEN_SPARSE_BLOCK_H
    410