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_SPARSE_BLOCK_H
11#define EIGEN_SPARSE_BLOCK_H
12
13namespace Eigen {
14
15template<typename XprType, int BlockRows, int BlockCols>
16class 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;
21public:
22    enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
23protected:
24    enum { OuterSize = IsRowMajor ? BlockRows : BlockCols };
25public:
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
78template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols>
79class 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;
85public:
86    enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
87    EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType)
88protected:
89    enum { OuterSize = IsRowMajor ? BlockRows : BlockCols };
90public:
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  */
253template<typename Derived>
254typename 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  */
260template<typename Derived>
261const 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  */
267template<typename Derived>
268Block<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  */
279template<typename Derived>
280const 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  */
291template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
292class 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;
297public:
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