DynamicSparseMatrix.h
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_DYNAMIC_SPARSEMATRIX_H
00026 #define EIGEN_DYNAMIC_SPARSEMATRIX_H
00027 
00028 namespace Eigen { 
00029 
00050 namespace internal {
00051 template<typename _Scalar, int _Options, typename _Index>
00052 struct traits<DynamicSparseMatrix<_Scalar, _Options, _Index> >
00053 {
00054   typedef _Scalar Scalar;
00055   typedef _Index Index;
00056   typedef Sparse StorageKind;
00057   typedef MatrixXpr XprKind;
00058   enum {
00059     RowsAtCompileTime = Dynamic,
00060     ColsAtCompileTime = Dynamic,
00061     MaxRowsAtCompileTime = Dynamic,
00062     MaxColsAtCompileTime = Dynamic,
00063     Flags = _Options | NestByRefBit | LvalueBit,
00064     CoeffReadCost = NumTraits<Scalar>::ReadCost,
00065     SupportedAccessPatterns = OuterRandomAccessPattern
00066   };
00067 };
00068 }
00069 
00070 template<typename _Scalar, int _Options, typename _Index>
00071  class  DynamicSparseMatrix
00072   : public SparseMatrixBase<DynamicSparseMatrix<_Scalar, _Options, _Index> >
00073 {
00074   public:
00075     EIGEN_SPARSE_PUBLIC_INTERFACE(DynamicSparseMatrix)
00076     // FIXME: why are these operator already alvailable ???
00077     // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, +=)
00078     // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=)
00079     typedef MappedSparseMatrix<Scalar,Flags> Map;
00080     using Base::IsRowMajor;
00081     using Base::operator=;
00082     enum {
00083       Options = _Options
00084     };
00085 
00086   protected:
00087 
00088     typedef DynamicSparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix;
00089 
00090     Index m_innerSize;
00091     std::vector<internal::CompressedStorage<Scalar,Index> > m_data;
00092 
00093   public:
00094 
00095     inline Index rows() const { return IsRowMajor ? outerSize() : m_innerSize; }
00096     inline Index cols() const { return IsRowMajor ? m_innerSize : outerSize(); }
00097     inline Index innerSize() const { return m_innerSize; }
00098     inline Index outerSize() const { return static_cast<Index>(m_data.size()); }
00099     inline Index innerNonZeros(Index j) const { return m_data[j].size(); }
00100 
00101     std::vector<internal::CompressedStorage<Scalar,Index> >& _data() { return m_data; }
00102     const std::vector<internal::CompressedStorage<Scalar,Index> >& _data() const { return m_data; }
00103 
00107     inline Scalar coeff(Index row, Index col) const
00108     {
00109       const Index outer = IsRowMajor ? row : col;
00110       const Index inner = IsRowMajor ? col : row;
00111       return m_data[outer].at(inner);
00112     }
00113 
00118     inline Scalar& coeffRef(Index row, Index col)
00119     {
00120       const Index outer = IsRowMajor ? row : col;
00121       const Index inner = IsRowMajor ? col : row;
00122       return m_data[outer].atWithInsertion(inner);
00123     }
00124 
00125     class InnerIterator;
00126     class ReverseInnerIterator;
00127 
00128     void setZero()
00129     {
00130       for (Index j=0; j<outerSize(); ++j)
00131         m_data[j].clear();
00132     }
00133 
00135     Index nonZeros() const
00136     {
00137       Index res = 0;
00138       for (Index j=0; j<outerSize(); ++j)
00139         res += static_cast<Index>(m_data[j].size());
00140       return res;
00141     }
00142 
00143 
00144 
00145     void reserve(Index reserveSize = 1000)
00146     {
00147       if (outerSize()>0)
00148       {
00149         Index reserveSizePerVector = (std::max)(reserveSize/outerSize(),Index(4));
00150         for (Index j=0; j<outerSize(); ++j)
00151         {
00152           m_data[j].reserve(reserveSizePerVector);
00153         }
00154       }
00155     }
00156 
00158     inline void startVec(Index /*outer*/) {}
00159 
00165     inline Scalar& insertBack(Index row, Index col)
00166     {
00167       return insertBackByOuterInner(IsRowMajor?row:col, IsRowMajor?col:row);
00168     }
00169 
00171     inline Scalar& insertBackByOuterInner(Index outer, Index inner)
00172     {
00173       eigen_assert(outer<Index(m_data.size()) && inner<m_innerSize && "out of range");
00174       eigen_assert(((m_data[outer].size()==0) || (m_data[outer].index(m_data[outer].size()-1)<inner))
00175                 && "wrong sorted insertion");
00176       m_data[outer].append(0, inner);
00177       return m_data[outer].value(m_data[outer].size()-1);
00178     }
00179 
00180     inline Scalar& insert(Index row, Index col)
00181     {
00182       const Index outer = IsRowMajor ? row : col;
00183       const Index inner = IsRowMajor ? col : row;
00184 
00185       Index startId = 0;
00186       Index id = static_cast<Index>(m_data[outer].size()) - 1;
00187       m_data[outer].resize(id+2,1);
00188 
00189       while ( (id >= startId) && (m_data[outer].index(id) > inner) )
00190       {
00191         m_data[outer].index(id+1) = m_data[outer].index(id);
00192         m_data[outer].value(id+1) = m_data[outer].value(id);
00193         --id;
00194       }
00195       m_data[outer].index(id+1) = inner;
00196       m_data[outer].value(id+1) = 0;
00197       return m_data[outer].value(id+1);
00198     }
00199 
00201     inline void finalize() {}
00202 
00204     void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
00205     {
00206       for (Index j=0; j<outerSize(); ++j)
00207         m_data[j].prune(reference,epsilon);
00208     }
00209 
00212     void resize(Index rows, Index cols)
00213     {
00214       const Index outerSize = IsRowMajor ? rows : cols;
00215       m_innerSize = IsRowMajor ? cols : rows;
00216       setZero();
00217       if (Index(m_data.size()) != outerSize)
00218       {
00219         m_data.resize(outerSize);
00220       }
00221     }
00222 
00223     void resizeAndKeepData(Index rows, Index cols)
00224     {
00225       const Index outerSize = IsRowMajor ? rows : cols;
00226       const Index innerSize = IsRowMajor ? cols : rows;
00227       if (m_innerSize>innerSize)
00228       {
00229         // remove all coefficients with innerCoord>=innerSize
00230         // TODO
00231         //std::cerr << "not implemented yet\n";
00232         exit(2);
00233       }
00234       if (m_data.size() != outerSize)
00235       {
00236         m_data.resize(outerSize);
00237       }
00238     }
00239 
00241     EIGEN_DEPRECATED inline DynamicSparseMatrix()
00242       : m_innerSize(0), m_data(0)
00243     {
00244       eigen_assert(innerSize()==0 && outerSize()==0);
00245     }
00246 
00248     EIGEN_DEPRECATED inline DynamicSparseMatrix(Index rows, Index cols)
00249       : m_innerSize(0)
00250     {
00251       resize(rows, cols);
00252     }
00253 
00255     template<typename OtherDerived>
00256     EIGEN_DEPRECATED explicit inline DynamicSparseMatrix(const SparseMatrixBase<OtherDerived>& other)
00257       : m_innerSize(0)
00258     {
00259     Base::operator=(other.derived());
00260     }
00261 
00262     inline DynamicSparseMatrix(const DynamicSparseMatrix& other)
00263       : Base(), m_innerSize(0)
00264     {
00265       *this = other.derived();
00266     }
00267 
00268     inline void swap(DynamicSparseMatrix& other)
00269     {
00270       //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
00271       std::swap(m_innerSize, other.m_innerSize);
00272       //std::swap(m_outerSize, other.m_outerSize);
00273       m_data.swap(other.m_data);
00274     }
00275 
00276     inline DynamicSparseMatrix& operator=(const DynamicSparseMatrix& other)
00277     {
00278       if (other.isRValue())
00279       {
00280         swap(other.const_cast_derived());
00281       }
00282       else
00283       {
00284         resize(other.rows(), other.cols());
00285         m_data = other.m_data;
00286       }
00287       return *this;
00288     }
00289 
00291     inline ~DynamicSparseMatrix() {}
00292 
00293   public:
00294 
00297     EIGEN_DEPRECATED void startFill(Index reserveSize = 1000)
00298     {
00299       setZero();
00300       reserve(reserveSize);
00301     }
00302 
00312     EIGEN_DEPRECATED Scalar& fill(Index row, Index col)
00313     {
00314       const Index outer = IsRowMajor ? row : col;
00315       const Index inner = IsRowMajor ? col : row;
00316       return insertBack(outer,inner);
00317     }
00318 
00324     EIGEN_DEPRECATED Scalar& fillrand(Index row, Index col)
00325     {
00326       return insert(row,col);
00327     }
00328 
00331     EIGEN_DEPRECATED void endFill() {}
00332     
00333 #   ifdef EIGEN_DYNAMICSPARSEMATRIX_PLUGIN
00334 #     include EIGEN_DYNAMICSPARSEMATRIX_PLUGIN
00335 #   endif
00336  };
00337 
00338 template<typename Scalar, int _Options, typename _Index>
00339 class DynamicSparseMatrix<Scalar,_Options,_Index>::InnerIterator : public SparseVector<Scalar,_Options,_Index>::InnerIterator
00340 {
00341     typedef typename SparseVector<Scalar,_Options,_Index>::InnerIterator Base;
00342   public:
00343     InnerIterator(const DynamicSparseMatrix& mat, Index outer)
00344       : Base(mat.m_data[outer]), m_outer(outer)
00345     {}
00346 
00347     inline Index row() const { return IsRowMajor ? m_outer : Base::index(); }
00348     inline Index col() const { return IsRowMajor ? Base::index() : m_outer; }
00349 
00350   protected:
00351     const Index m_outer;
00352 };
00353 
00354 template<typename Scalar, int _Options, typename _Index>
00355 class DynamicSparseMatrix<Scalar,_Options,_Index>::ReverseInnerIterator : public SparseVector<Scalar,_Options,_Index>::ReverseInnerIterator
00356 {
00357     typedef typename SparseVector<Scalar,_Options,_Index>::ReverseInnerIterator Base;
00358   public:
00359     ReverseInnerIterator(const DynamicSparseMatrix& mat, Index outer)
00360       : Base(mat.m_data[outer]), m_outer(outer)
00361     {}
00362 
00363     inline Index row() const { return IsRowMajor ? m_outer : Base::index(); }
00364     inline Index col() const { return IsRowMajor ? Base::index() : m_outer; }
00365 
00366   protected:
00367     const Index m_outer;
00368 };
00369 
00370 } // end namespace Eigen
00371 
00372 #endif // EIGEN_DYNAMIC_SPARSEMATRIX_H