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 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_SPARSE_BLOCK_H
00011 #define EIGEN_SPARSE_BLOCK_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 template<typename MatrixType, int Size>
00017 struct traits<SparseInnerVectorSet<MatrixType, Size> >
00018 {
00019   typedef typename traits<MatrixType>::Scalar Scalar;
00020   typedef typename traits<MatrixType>::Index Index;
00021   typedef typename traits<MatrixType>::StorageKind StorageKind;
00022   typedef MatrixXpr XprKind;
00023   enum {
00024     IsRowMajor = (int(MatrixType::Flags)&RowMajorBit)==RowMajorBit,
00025     Flags = MatrixType::Flags,
00026     RowsAtCompileTime = IsRowMajor ? Size : MatrixType::RowsAtCompileTime,
00027     ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : Size,
00028     MaxRowsAtCompileTime = RowsAtCompileTime,
00029     MaxColsAtCompileTime = ColsAtCompileTime,
00030     CoeffReadCost = MatrixType::CoeffReadCost
00031   };
00032 };
00033 } // end namespace internal
00034 
00035 template<typename MatrixType, int Size>
00036 class SparseInnerVectorSet : internal::no_assignment_operator,
00037   public SparseMatrixBase<SparseInnerVectorSet<MatrixType, Size> >
00038 {
00039   public:
00040 
00041     enum { IsRowMajor = internal::traits<SparseInnerVectorSet>::IsRowMajor };
00042 
00043     EIGEN_SPARSE_PUBLIC_INTERFACE(SparseInnerVectorSet)
00044     class InnerIterator: public MatrixType::InnerIterator
00045     {
00046       public:
00047         inline InnerIterator(const SparseInnerVectorSet& xpr, Index outer)
00048           : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
00049         {}
00050         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
00051         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
00052       protected:
00053         Index m_outer;
00054     };
00055     class ReverseInnerIterator: public MatrixType::ReverseInnerIterator
00056     {
00057       public:
00058         inline ReverseInnerIterator(const SparseInnerVectorSet& xpr, Index outer)
00059           : MatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
00060         {}
00061         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
00062         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
00063       protected:
00064         Index m_outer;
00065     };
00066 
00067     inline SparseInnerVectorSet(const MatrixType& matrix, Index outerStart, Index outerSize)
00068       : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
00069     {
00070       eigen_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
00071     }
00072 
00073     inline SparseInnerVectorSet(const MatrixType& matrix, Index outer)
00074       : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
00075     {
00076       eigen_assert(Size!=Dynamic);
00077       eigen_assert( (outer>=0) && (outer<matrix.outerSize()) );
00078     }
00079 
00080 //     template<typename OtherDerived>
00081 //     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
00082 //     {
00083 //       return *this;
00084 //     }
00085 
00086 //     template<typename Sparse>
00087 //     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
00088 //     {
00089 //       return *this;
00090 //     }
00091 
00092     EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
00093     EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
00094 
00095   protected:
00096 
00097     const typename MatrixType::Nested m_matrix;
00098     Index m_outerStart;
00099     const internal::variable_if_dynamic<Index, Size> m_outerSize;
00100 };
00101 
00102 
00103 /***************************************************************************
00104 * specialisation for SparseMatrix
00105 ***************************************************************************/
00106 
00107 template<typename _Scalar, int _Options, typename _Index, int Size>
00108 class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options, _Index>, Size>
00109   : public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options, _Index>, Size> >
00110 {
00111     typedef SparseMatrix<_Scalar, _Options, _Index> MatrixType;
00112   public:
00113 
00114     enum { IsRowMajor = internal::traits<SparseInnerVectorSet>::IsRowMajor };
00115 
00116     EIGEN_SPARSE_PUBLIC_INTERFACE(SparseInnerVectorSet)
00117     class InnerIterator: public MatrixType::InnerIterator
00118     {
00119       public:
00120         inline InnerIterator(const SparseInnerVectorSet& xpr, Index outer)
00121           : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
00122         {}
00123         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
00124         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
00125       protected:
00126         Index m_outer;
00127     };
00128     class ReverseInnerIterator: public MatrixType::ReverseInnerIterator
00129     {
00130       public:
00131         inline ReverseInnerIterator(const SparseInnerVectorSet& xpr, Index outer)
00132           : MatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
00133         {}
00134         inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
00135         inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
00136       protected:
00137         Index m_outer;
00138     };
00139 
00140     inline SparseInnerVectorSet(const MatrixType& matrix, Index outerStart, Index outerSize)
00141       : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
00142     {
00143       eigen_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
00144     }
00145 
00146     inline SparseInnerVectorSet(const MatrixType& matrix, Index outer)
00147       : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
00148     {
00149       eigen_assert(Size==1);
00150       eigen_assert( (outer>=0) && (outer<matrix.outerSize()) );
00151     }
00152 
00153     template<typename OtherDerived>
00154     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
00155     {
00156       typedef typename internal::remove_all<typename MatrixType::Nested>::type _NestedMatrixType;
00157       _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);;
00158       // This assignement is slow if this vector set is not empty
00159       // and/or it is not at the end of the nonzeros of the underlying matrix.
00160 
00161       // 1 - eval to a temporary to avoid transposition and/or aliasing issues
00162       SparseMatrix<Scalar, IsRowMajor ? RowMajor : ColMajor, Index> tmp(other);
00163 
00164       // 2 - let's check whether there is enough allocated memory
00165       Index nnz           = tmp.nonZeros();
00166       Index nnz_previous  = nonZeros();
00167       Index free_size     = Index(matrix.data().allocatedSize()) + nnz_previous;
00168       Index nnz_head      = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart];
00169       Index tail          = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()];
00170       Index nnz_tail      = matrix.nonZeros() - tail;
00171 
00172       if(nnz>free_size)
00173       {
00174         // realloc manually to reduce copies
00175         typename MatrixType::Storage newdata(m_matrix.nonZeros() - nnz_previous + nnz);
00176 
00177         std::memcpy(&newdata.value(0), &m_matrix.data().value(0), nnz_head*sizeof(Scalar));
00178         std::memcpy(&newdata.index(0), &m_matrix.data().index(0), nnz_head*sizeof(Index));
00179 
00180         std::memcpy(&newdata.value(nnz_head), &tmp.data().value(0), nnz*sizeof(Scalar));
00181         std::memcpy(&newdata.index(nnz_head), &tmp.data().index(0), nnz*sizeof(Index));
00182 
00183         std::memcpy(&newdata.value(nnz_head+nnz), &matrix.data().value(tail), nnz_tail*sizeof(Scalar));
00184         std::memcpy(&newdata.index(nnz_head+nnz), &matrix.data().index(tail), nnz_tail*sizeof(Index));
00185 
00186         matrix.data().swap(newdata);
00187       }
00188       else
00189       {
00190         // no need to realloc, simply copy the tail at its respective position and insert tmp
00191         matrix.data().resize(nnz_head + nnz + nnz_tail);
00192 
00193         if(nnz<nnz_previous)
00194         {
00195           std::memcpy(&matrix.data().value(nnz_head+nnz), &matrix.data().value(tail), nnz_tail*sizeof(Scalar));
00196           std::memcpy(&matrix.data().index(nnz_head+nnz), &matrix.data().index(tail), nnz_tail*sizeof(Index));
00197         }
00198         else
00199         {
00200           for(Index i=nnz_tail-1; i>=0; --i)
00201           {
00202             matrix.data().value(nnz_head+nnz+i) = matrix.data().value(tail+i);
00203             matrix.data().index(nnz_head+nnz+i) = matrix.data().index(tail+i);
00204           }
00205         }
00206 
00207         std::memcpy(&matrix.data().value(nnz_head), &tmp.data().value(0), nnz*sizeof(Scalar));
00208         std::memcpy(&matrix.data().index(nnz_head), &tmp.data().index(0), nnz*sizeof(Index));
00209       }
00210 
00211       // update outer index pointers
00212       Index p = nnz_head;
00213       for(Index k=0; k<m_outerSize.value(); ++k)
00214       {
00215         matrix.outerIndexPtr()[m_outerStart+k] = p;
00216         p += tmp.innerVector(k).nonZeros();
00217       }
00218       std::ptrdiff_t offset = nnz - nnz_previous;
00219       for(Index k = m_outerStart + m_outerSize.value(); k<=matrix.outerSize(); ++k)
00220       {
00221         matrix.outerIndexPtr()[k] += offset;
00222       }
00223 
00224       return *this;
00225     }
00226 
00227     inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other)
00228     {
00229       return operator=<SparseInnerVectorSet>(other);
00230     }
00231 
00232     inline const Scalar* valuePtr() const
00233     { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
00234     inline Scalar* valuePtr()
00235     { return m_matrix.const_cast_derived().valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
00236 
00237     inline const Index* innerIndexPtr() const
00238     { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
00239     inline Index* innerIndexPtr()
00240     { return m_matrix.const_cast_derived().innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
00241 
00242     inline const Index* outerIndexPtr() const
00243     { return m_matrix.outerIndexPtr() + m_outerStart; }
00244     inline Index* outerIndexPtr()
00245     { return m_matrix.const_cast_derived().outerIndexPtr() + m_outerStart; }
00246 
00247     Index nonZeros() const
00248     {
00249       if(m_matrix.isCompressed())
00250         return  std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()])
00251               - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]);
00252       else if(m_outerSize.value()==0)
00253         return 0;
00254       else
00255         return Map<const Matrix<Index,Size,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum();
00256     }
00257 
00258     const Scalar& lastCoeff() const
00259     {
00260       EIGEN_STATIC_ASSERT_VECTOR_ONLY(SparseInnerVectorSet);
00261       eigen_assert(nonZeros()>0);
00262       if(m_matrix.isCompressed())
00263         return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1];
00264       else
00265         return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1];
00266     }
00267 
00268 //     template<typename Sparse>
00269 //     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
00270 //     {
00271 //       return *this;
00272 //     }
00273 
00274     EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
00275     EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
00276 
00277   protected:
00278 
00279     typename MatrixType::Nested m_matrix;
00280     Index m_outerStart;
00281     const internal::variable_if_dynamic<Index, Size> m_outerSize;
00282 
00283 };
00284 
00285 //----------
00286 
00288 template<typename Derived>
00289 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(Index i)
00290 {
00291   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
00292   return innerVector(i);
00293 }
00294 
00297 template<typename Derived>
00298 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(Index i) const
00299 {
00300   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
00301   return innerVector(i);
00302 }
00303 
00305 template<typename Derived>
00306 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(Index i)
00307 {
00308   EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00309   return innerVector(i);
00310 }
00311 
00314 template<typename Derived>
00315 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(Index i) const
00316 {
00317   EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00318   return innerVector(i);
00319 }
00320 
00324 template<typename Derived>
00325 SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(Index outer)
00326 { return SparseInnerVectorSet<Derived,1>(derived(), outer); }
00327 
00331 template<typename Derived>
00332 const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(Index outer) const
00333 { return SparseInnerVectorSet<Derived,1>(derived(), outer); }
00334 
00336 template<typename Derived>
00337 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleRows(Index start, Index size)
00338 {
00339   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
00340   return innerVectors(start, size);
00341 }
00342 
00345 template<typename Derived>
00346 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleRows(Index start, Index size) const
00347 {
00348   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
00349   return innerVectors(start, size);
00350 }
00351 
00353 template<typename Derived>
00354 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleCols(Index start, Index size)
00355 {
00356   EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00357   return innerVectors(start, size);
00358 }
00359 
00362 template<typename Derived>
00363 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::middleCols(Index start, Index size) const
00364 {
00365   EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00366   return innerVectors(start, size);
00367 }
00368 
00369 
00370 
00374 template<typename Derived>
00375 SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize)
00376 { return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
00377 
00381 template<typename Derived>
00382 const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize) const
00383 { return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
00384 
00385 } // end namespace Eigen
00386 
00387 #endif // EIGEN_SPARSE_BLOCK_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:26:16