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