SparseBlock.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_SPARSE_BLOCK_H
00026 #define EIGEN_SPARSE_BLOCK_H
00027 
00028 namespace Eigen { 
00029 
00030 namespace internal {
00031 template<typename MatrixType, int Size>
00032 struct traits<SparseInnerVectorSet<MatrixType, Size> >
00033 {
00034   typedef typename traits<MatrixType>::Scalar Scalar;
00035   typedef typename traits<MatrixType>::Index Index;
00036   typedef typename traits<MatrixType>::StorageKind StorageKind;
00037   typedef MatrixXpr XprKind;
00038   enum {
00039     IsRowMajor = (int(MatrixType::Flags)&RowMajorBit)==RowMajorBit,
00040     Flags = MatrixType::Flags,
00041     RowsAtCompileTime = IsRowMajor ? Size : MatrixType::RowsAtCompileTime,
00042     ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : Size,
00043     MaxRowsAtCompileTime = RowsAtCompileTime,
00044     MaxColsAtCompileTime = ColsAtCompileTime,
00045     CoeffReadCost = MatrixType::CoeffReadCost
00046   };
00047 };
00048 } // end namespace internal
00049 
00050 template<typename MatrixType, int Size>
00051 class SparseInnerVectorSet : internal::no_assignment_operator,
00052   public SparseMatrixBase<SparseInnerVectorSet<MatrixType, Size> >
00053 {
00054   public:
00055 
00056     enum { IsRowMajor = internal::traits<SparseInnerVectorSet>::IsRowMajor };
00057 
00058     EIGEN_SPARSE_PUBLIC_INTERFACE(SparseInnerVectorSet)
00059     class InnerIterator: public MatrixType::InnerIterator
00060     {
00061       public:
00062         inline InnerIterator(const SparseInnerVectorSet& xpr, Index outer)
00063           : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
00064         {}
00065         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
00066         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
00067       protected:
00068         Index m_outer;
00069     };
00070     class ReverseInnerIterator: public MatrixType::ReverseInnerIterator
00071     {
00072       public:
00073         inline ReverseInnerIterator(const SparseInnerVectorSet& xpr, Index outer)
00074           : MatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
00075         {}
00076         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
00077         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
00078       protected:
00079         Index m_outer;
00080     };
00081 
00082     inline SparseInnerVectorSet(const MatrixType& matrix, Index outerStart, Index outerSize)
00083       : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
00084     {
00085       eigen_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
00086     }
00087 
00088     inline SparseInnerVectorSet(const MatrixType& matrix, Index outer)
00089       : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
00090     {
00091       eigen_assert(Size!=Dynamic);
00092       eigen_assert( (outer>=0) && (outer<matrix.outerSize()) );
00093     }
00094 
00095 //     template<typename OtherDerived>
00096 //     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
00097 //     {
00098 //       return *this;
00099 //     }
00100 
00101 //     template<typename Sparse>
00102 //     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
00103 //     {
00104 //       return *this;
00105 //     }
00106 
00107     EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
00108     EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
00109 
00110   protected:
00111 
00112     const typename MatrixType::Nested m_matrix;
00113     Index m_outerStart;
00114     const internal::variable_if_dynamic<Index, Size> m_outerSize;
00115 };
00116 
00117 
00118 /***************************************************************************
00119 * specialisation for SparseMatrix
00120 ***************************************************************************/
00121 
00122 template<typename _Scalar, int _Options, typename _Index, int Size>
00123 class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options, _Index>, Size>
00124   : public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options, _Index>, Size> >
00125 {
00126     typedef SparseMatrix<_Scalar, _Options, _Index> MatrixType;
00127   public:
00128 
00129     enum { IsRowMajor = internal::traits<SparseInnerVectorSet>::IsRowMajor };
00130 
00131     EIGEN_SPARSE_PUBLIC_INTERFACE(SparseInnerVectorSet)
00132     class InnerIterator: public MatrixType::InnerIterator
00133     {
00134       public:
00135         inline InnerIterator(const SparseInnerVectorSet& xpr, Index outer)
00136           : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
00137         {}
00138         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
00139         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
00140       protected:
00141         Index m_outer;
00142     };
00143     class ReverseInnerIterator: public MatrixType::ReverseInnerIterator
00144     {
00145       public:
00146         inline ReverseInnerIterator(const SparseInnerVectorSet& xpr, Index outer)
00147           : MatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
00148         {}
00149         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
00150         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
00151       protected:
00152         Index m_outer;
00153     };
00154 
00155     inline SparseInnerVectorSet(const MatrixType& matrix, Index outerStart, Index outerSize)
00156       : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
00157     {
00158       eigen_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
00159     }
00160 
00161     inline SparseInnerVectorSet(const MatrixType& matrix, Index outer)
00162       : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
00163     {
00164       eigen_assert(Size==1);
00165       eigen_assert( (outer>=0) && (outer<matrix.outerSize()) );
00166     }
00167 
00168     template<typename OtherDerived>
00169     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
00170     {
00171       typedef typename internal::remove_all<typename MatrixType::Nested>::type _NestedMatrixType;
00172       _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);;
00173       // This assignement is slow if this vector set is not empty
00174       // and/or it is not at the end of the nonzeros of the underlying matrix.
00175 
00176       // 1 - eval to a temporary to avoid transposition and/or aliasing issues
00177       SparseMatrix<Scalar, IsRowMajor ? RowMajor : ColMajor, Index> tmp(other);
00178 
00179       // 2 - let's check whether there is enough allocated memory
00180       Index nnz           = tmp.nonZeros();
00181       Index nnz_previous  = nonZeros();
00182       Index free_size     = Index(matrix.data().allocatedSize()) + nnz_previous;
00183       Index nnz_head      = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart];
00184       Index tail          = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()];
00185       Index nnz_tail      = matrix.nonZeros() - tail;
00186 
00187       if(nnz>free_size)
00188       {
00189         // realloc manually to reduce copies
00190         typename MatrixType::Storage newdata(m_matrix.nonZeros() - nnz_previous + nnz);
00191 
00192         std::memcpy(&newdata.value(0), &m_matrix.data().value(0), nnz_head*sizeof(Scalar));
00193         std::memcpy(&newdata.index(0), &m_matrix.data().index(0), nnz_head*sizeof(Index));
00194 
00195         std::memcpy(&newdata.value(nnz_head), &tmp.data().value(0), nnz*sizeof(Scalar));
00196         std::memcpy(&newdata.index(nnz_head), &tmp.data().index(0), nnz*sizeof(Index));
00197 
00198         std::memcpy(&newdata.value(nnz_head+nnz), &matrix.data().value(tail), nnz_tail*sizeof(Scalar));
00199         std::memcpy(&newdata.index(nnz_head+nnz), &matrix.data().index(tail), nnz_tail*sizeof(Index));
00200 
00201         matrix.data().swap(newdata);
00202       }
00203       else
00204       {
00205         // no need to realloc, simply copy the tail at its respective position and insert tmp
00206         matrix.data().resize(nnz_head + nnz + nnz_tail);
00207 
00208         if(nnz<nnz_previous)
00209         {
00210           std::memcpy(&matrix.data().value(nnz_head+nnz), &matrix.data().value(tail), nnz_tail*sizeof(Scalar));
00211           std::memcpy(&matrix.data().index(nnz_head+nnz), &matrix.data().index(tail), nnz_tail*sizeof(Index));
00212         }
00213         else
00214         {
00215           for(Index i=nnz_tail-1; i>=0; --i)
00216           {
00217             matrix.data().value(nnz_head+nnz+i) = matrix.data().value(tail+i);
00218             matrix.data().index(nnz_head+nnz+i) = matrix.data().index(tail+i);
00219           }
00220         }
00221 
00222         std::memcpy(&matrix.data().value(nnz_head), &tmp.data().value(0), nnz*sizeof(Scalar));
00223         std::memcpy(&matrix.data().index(nnz_head), &tmp.data().index(0), nnz*sizeof(Index));
00224       }
00225 
00226       // update outer index pointers
00227       Index p = nnz_head;
00228       for(Index k=0; k<m_outerSize.value(); ++k)
00229       {
00230         matrix.outerIndexPtr()[m_outerStart+k] = p;
00231         p += tmp.innerVector(k).nonZeros();
00232       }
00233       std::ptrdiff_t offset = nnz - nnz_previous;
00234       for(Index k = m_outerStart + m_outerSize.value(); k<=matrix.outerSize(); ++k)
00235       {
00236         matrix.outerIndexPtr()[k] += offset;
00237       }
00238 
00239       return *this;
00240     }
00241 
00242     inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other)
00243     {
00244       return operator=<SparseInnerVectorSet>(other);
00245     }
00246 
00247     inline const Scalar* valuePtr() const
00248     { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
00249     inline Scalar* valuePtr()
00250     { return m_matrix.const_cast_derived().valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
00251 
00252     inline const Index* innerIndexPtr() const
00253     { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
00254     inline Index* innerIndexPtr()
00255     { return m_matrix.const_cast_derived().innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
00256 
00257     inline const Index* outerIndexPtr() const
00258     { return m_matrix.outerIndexPtr() + m_outerStart; }
00259     inline Index* outerIndexPtr()
00260     { return m_matrix.const_cast_derived().outerIndexPtr() + m_outerStart; }
00261 
00262     Index nonZeros() const
00263     {
00264       if(m_matrix.isCompressed())
00265         return  std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()])
00266               - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]);
00267       else if(m_outerSize.value()==0)
00268         return 0;
00269       else
00270         return Map<const Matrix<Index,Size,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum();
00271     }
00272 
00273     const Scalar& lastCoeff() const
00274     {
00275       EIGEN_STATIC_ASSERT_VECTOR_ONLY(SparseInnerVectorSet);
00276       eigen_assert(nonZeros()>0);
00277       if(m_matrix.isCompressed())
00278         return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1];
00279       else
00280         return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1];
00281     }
00282 
00283 //     template<typename Sparse>
00284 //     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
00285 //     {
00286 //       return *this;
00287 //     }
00288 
00289     EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
00290     EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
00291 
00292   protected:
00293 
00294     typename MatrixType::Nested m_matrix;
00295     Index m_outerStart;
00296     const internal::variable_if_dynamic<Index, Size> m_outerSize;
00297 
00298 };
00299 
00300 //----------
00301 
00303 template<typename Derived>
00304 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(Index i)
00305 {
00306   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
00307   return innerVector(i);
00308 }
00309 
00312 template<typename Derived>
00313 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(Index i) const
00314 {
00315   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
00316   return innerVector(i);
00317 }
00318 
00320 template<typename Derived>
00321 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(Index i)
00322 {
00323   EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00324   return innerVector(i);
00325 }
00326 
00329 template<typename Derived>
00330 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(Index i) const
00331 {
00332   EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00333   return innerVector(i);
00334 }
00335 
00339 template<typename Derived>
00340 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(Index outer)
00341 { return SparseInnerVectorSet<Derived,1>(derived(), outer); }
00342 
00346 template<typename Derived>
00347 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(Index outer) const
00348 { return SparseInnerVectorSet<Derived,1>(derived(), outer); }
00349 
00351 template<typename Derived>
00352 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleRows(Index start, Index size)
00353 {
00354   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
00355   return innerVectors(start, size);
00356 }
00357 
00360 template<typename Derived>
00361 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleRows(Index start, Index size) const
00362 {
00363   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
00364   return innerVectors(start, size);
00365 }
00366 
00368 template<typename Derived>
00369 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleCols(Index start, Index size)
00370 {
00371   EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00372   return innerVectors(start, size);
00373 }
00374 
00377 template<typename Derived>
00378 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleCols(Index start, Index size) const
00379 {
00380   EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00381   return innerVectors(start, size);
00382 }
00383 
00384 
00385 
00389 template<typename Derived>
00390 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize)
00391 { return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
00392 
00396 template<typename Derived>
00397 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize) const
00398 { return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
00399 
00400 } // end namespace Eigen
00401 
00402 #endif // EIGEN_SPARSE_BLOCK_H