diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-05-20 00:05:28 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-05-20 00:05:28 +0200 |
commit | 6ecd02d7ec85f07e02559cb311d4dd07e844a72d (patch) | |
tree | 3703e87d24166af9d2e74a8522c92c3393b6557b /Eigen/src | |
parent | 56aad5aafb84b87c9cdf8a237a3d7fdf632b2021 (diff) | |
parent | a697cb409440c3d01eda44f84abe2f1d924270ce (diff) |
merge
Diffstat (limited to 'Eigen/src')
-rw-r--r-- | Eigen/src/LU/PartialLU.h | 4 | ||||
-rw-r--r-- | Eigen/src/Sparse/DynamicSparseMatrix.h | 30 | ||||
-rw-r--r-- | Eigen/src/Sparse/SparseBlock.h | 27 | ||||
-rw-r--r-- | Eigen/src/Sparse/SparseCwiseBinaryOp.h | 4 | ||||
-rw-r--r-- | Eigen/src/Sparse/SparseMatrix.h | 32 | ||||
-rw-r--r-- | Eigen/src/Sparse/SparseMatrixBase.h | 35 | ||||
-rw-r--r-- | Eigen/src/Sparse/SparseRedux.h | 16 | ||||
-rw-r--r-- | Eigen/src/Sparse/SparseTranspose.h | 10 | ||||
-rw-r--r-- | Eigen/src/Sparse/SparseVector.h | 27 |
9 files changed, 125 insertions, 60 deletions
diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h index 7fdbeac38..3a6e9f286 100644 --- a/Eigen/src/LU/PartialLU.h +++ b/Eigen/src/LU/PartialLU.h @@ -225,8 +225,8 @@ void PartialLU<MatrixType>::solve( /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. * So we proceed as follows: * Step 1: compute c = Pb. - * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible. - * Step 3: replace c by the solution x to Ux = c. Check if a solution really exists. + * Step 2: replace c by the solution x to Lx = c. + * Step 3: replace c by the solution x to Ux = c. */ const int size = m_lu.rows(); diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h index 2927cd583..ec730dd3f 100644 --- a/Eigen/src/Sparse/DynamicSparseMatrix.h +++ b/Eigen/src/Sparse/DynamicSparseMatrix.h @@ -35,7 +35,7 @@ * random read/write accesses in log(rho*outer_size) where \c rho is the probability that a coefficient is * nonzero and outer_size is the number of columns if the matrix is column-major and the number of rows * otherwise. - * + * * Internally, the data are stored as a std::vector of compressed vector. The performances of random writes might * decrease as the number of nonzeros per inner-vector increase. In practice, we observed very good performance * till about 100 nonzeros/vector, and the performance remains relatively good till 500 nonzeros/vectors. @@ -67,23 +67,23 @@ class DynamicSparseMatrix // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, +=) // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=) typedef MappedSparseMatrix<Scalar,Flags> Map; + using Base::IsRowMajor; protected: - enum { IsRowMajor = Base::IsRowMajor }; typedef DynamicSparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix; int m_innerSize; std::vector<CompressedStorage<Scalar> > m_data; public: - + inline int rows() const { return IsRowMajor ? outerSize() : m_innerSize; } inline int cols() const { return IsRowMajor ? m_innerSize : outerSize(); } inline int innerSize() const { return m_innerSize; } inline int outerSize() const { return m_data.size(); } inline int innerNonZeros(int j) const { return m_data[j].size(); } - + std::vector<CompressedStorage<Scalar> >& _data() { return m_data; } const std::vector<CompressedStorage<Scalar> >& _data() const { return m_data; } @@ -132,7 +132,7 @@ class DynamicSparseMatrix setZero(); reserve(reserveSize); } - + void reserve(int reserveSize = 1000) { if (outerSize()>0) @@ -144,13 +144,13 @@ class DynamicSparseMatrix } } } - + inline void startVec(int /*outer*/) {} - + inline Scalar& insertBack(int outer, int inner) { ei_assert(outer<int(m_data.size()) && inner<m_innerSize && "out of range"); - ei_assert(((m_data[outer].size()==0) || (m_data[outer].index(m_data[outer].size()-1)<inner)) + ei_assert(((m_data[outer].size()==0) || (m_data[outer].index(m_data[outer].size()-1)<inner)) && "wrong sorted insertion"); m_data[outer].append(0, inner); return m_data[outer].value(m_data[outer].size()-1); @@ -162,7 +162,7 @@ class DynamicSparseMatrix * 2 - this the coefficient with greater inner coordinate for the given outer coordinate. * In other words, assuming \c *this is column-major, then there must not exists any nonzero coefficient of coordinates * \c i \c x \a col such that \c i >= \a row. Otherwise the matrix is invalid. - * + * * \see fillrand(), coeffRef() */ EIGEN_DEPRECATED Scalar& fill(int row, int col) @@ -181,12 +181,12 @@ class DynamicSparseMatrix { return insert(row,col); } - + inline Scalar& insert(int row, int col) { const int outer = IsRowMajor ? row : col; const int inner = IsRowMajor ? col : row; - + int startId = 0; int id = m_data[outer].size() - 1; m_data[outer].resize(id+2,1); @@ -205,9 +205,9 @@ class DynamicSparseMatrix /** \deprecated use finalize() * Does nothing. Provided for compatibility with SparseMatrix. */ EIGEN_DEPRECATED void endFill() {} - + inline void finalize() {} - + void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) { for (int j=0; j<outerSize(); ++j) @@ -226,7 +226,7 @@ class DynamicSparseMatrix m_data.resize(outerSize); } } - + void resizeAndKeepData(int rows, int cols) { const int outerSize = IsRowMajor ? rows : cols; @@ -309,7 +309,7 @@ class DynamicSparseMatrix<Scalar,_Flags>::InnerIterator : public SparseVector<Sc InnerIterator(const DynamicSparseMatrix& mat, int outer) : Base(mat.m_data[outer]), m_outer(outer) {} - + inline int row() const { return IsRowMajor ? m_outer : Base::index(); } inline int col() const { return IsRowMajor ? Base::index() : m_outer; } diff --git a/Eigen/src/Sparse/SparseBlock.h b/Eigen/src/Sparse/SparseBlock.h index 72589d0c0..655ba6d93 100644 --- a/Eigen/src/Sparse/SparseBlock.h +++ b/Eigen/src/Sparse/SparseBlock.h @@ -43,16 +43,21 @@ template<typename MatrixType, int Size> class SparseInnerVectorSet : ei_no_assignment_operator, public SparseMatrixBase<SparseInnerVectorSet<MatrixType, Size> > { - enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor }; public: + enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor }; + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet) class InnerIterator: public MatrixType::InnerIterator { public: inline InnerIterator(const SparseInnerVectorSet& xpr, int outer) - : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer) + : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) {} + inline int row() const { return IsRowMajor ? m_outer : this->index(); } + inline int col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + int m_outer; }; inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize) @@ -100,16 +105,21 @@ class SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size> : public SparseMatrixBase<SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size> > { typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType; - enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor }; public: + enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor }; + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet) class InnerIterator: public MatrixType::InnerIterator { public: inline InnerIterator(const SparseInnerVectorSet& xpr, int outer) - : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer) + : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) {} + inline int row() const { return IsRowMajor ? m_outer : this->index(); } + inline int col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + int m_outer; }; inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize) @@ -193,16 +203,21 @@ class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> : public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> > { typedef SparseMatrix<_Scalar, _Options> MatrixType; - enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor }; public: + enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor }; + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet) class InnerIterator: public MatrixType::InnerIterator { public: inline InnerIterator(const SparseInnerVectorSet& xpr, int outer) - : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer) + : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) {} + inline int row() const { return IsRowMajor ? m_outer : this->index(); } + inline int col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + int m_outer; }; inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize) diff --git a/Eigen/src/Sparse/SparseCwiseBinaryOp.h b/Eigen/src/Sparse/SparseCwiseBinaryOp.h index d19970efc..96fd1d36a 100644 --- a/Eigen/src/Sparse/SparseCwiseBinaryOp.h +++ b/Eigen/src/Sparse/SparseCwiseBinaryOp.h @@ -186,8 +186,8 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Deri EIGEN_STRONG_INLINE Scalar value() const { return m_value; } EIGEN_STRONG_INLINE int index() const { return m_id; } - EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); } - EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); } + EIGEN_STRONG_INLINE int row() const { return Lhs::IsRowMajor ? m_lhsIter.row() : index(); } + EIGEN_STRONG_INLINE int col() const { return Lhs::IsRowMajor ? index() : m_lhsIter.col(); } EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; } diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 8bcf3c8b9..65cee69cf 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -25,17 +25,24 @@ #ifndef EIGEN_SPARSEMATRIX_H #define EIGEN_SPARSEMATRIX_H -/** \class SparseMatrix +/** \ingroup Sparse_Module * - * \brief Sparse matrix + * \class SparseMatrix + * + * \brief The main sparse matrix class + * + * This class implements a sparse matrix using the very common compressed row/column storage + * scheme. * * \param _Scalar the scalar type, i.e. the type of the coefficients + * \param _Options Union of bit flags controlling the storage scheme. Currently the only possibility + * is RowMajor. The default is 0 which means column-major. * * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. * */ -template<typename _Scalar, int _Flags> -struct ei_traits<SparseMatrix<_Scalar, _Flags> > +template<typename _Scalar, int _Options> +struct ei_traits<SparseMatrix<_Scalar, _Options> > { typedef _Scalar Scalar; enum { @@ -43,17 +50,15 @@ struct ei_traits<SparseMatrix<_Scalar, _Flags> > ColsAtCompileTime = Dynamic, MaxRowsAtCompileTime = Dynamic, MaxColsAtCompileTime = Dynamic, - Flags = SparseBit | _Flags, + Flags = SparseBit | _Options, CoeffReadCost = NumTraits<Scalar>::ReadCost, SupportedAccessPatterns = InnerRandomAccessPattern }; }; - - -template<typename _Scalar, int _Flags> +template<typename _Scalar, int _Options> class SparseMatrix - : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> > + : public SparseMatrixBase<SparseMatrix<_Scalar, _Options> > { public: EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseMatrix) @@ -64,10 +69,10 @@ class SparseMatrix // EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, /=) typedef MappedSparseMatrix<Scalar,Flags> Map; + using Base::IsRowMajor; protected: - enum { IsRowMajor = Base::IsRowMajor }; typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix; int m_outerSize; @@ -508,10 +513,13 @@ class SparseMatrix { delete[] m_outerIndex; } + + /** Overloaded for performance */ + Scalar sum() const; }; -template<typename Scalar, int _Flags> -class SparseMatrix<Scalar,_Flags>::InnerIterator +template<typename Scalar, int _Options> +class SparseMatrix<Scalar,_Options>::InnerIterator { public: InnerIterator(const SparseMatrix& mat, int outer) diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 948424f90..34779d9e6 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -25,6 +25,17 @@ #ifndef EIGEN_SPARSEMATRIXBASE_H #define EIGEN_SPARSEMATRIXBASE_H +/** \ingroup Sparse_Module + * + * \class SparseMatrixBase + * + * \brief Base class of any sparse matrices or sparse expressions + * + * \param Derived + * + * + * + */ template<typename Derived> class SparseMatrixBase { public: @@ -69,7 +80,7 @@ template<typename Derived> class SparseMatrixBase /**< This stores expression \ref flags flags which may or may not be inherited by new expressions * constructed from this one. See the \ref flags "list of flags". */ - + CoeffReadCost = ei_traits<Derived>::CoeffReadCost, /**< This is a rough measure of how expensive it is to read one coefficient from * this expression. @@ -156,9 +167,9 @@ template<typename Derived> class SparseMatrixBase ei_assert(( ((ei_traits<Derived>::SupportedAccessPatterns&OuterRandomAccessPattern)==OuterRandomAccessPattern) || (!((Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit)))) && "the transpose operation is supposed to be handled in SparseMatrix::operator="); - + enum { Flip = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit) }; - + const int outerSize = other.outerSize(); //typedef typename ei_meta_if<transpose, LinkedVectorMatrix<Scalar,Flags&RowMajorBit>, Derived>::ret TempType; // thanks to shallow copies, we always eval to a tempary @@ -293,13 +304,13 @@ template<typename Derived> class SparseMatrixBase template<typename OtherDerived> const typename SparseProductReturnType<Derived,OtherDerived>::Type operator*(const SparseMatrixBase<OtherDerived> &other) const; - + // dense * sparse (return a dense object) - template<typename OtherDerived> friend + template<typename OtherDerived> friend const typename SparseProductReturnType<OtherDerived,Derived>::Type operator*(const MatrixBase<OtherDerived>& lhs, const Derived& rhs) { return typename SparseProductReturnType<OtherDerived,Derived>::Type(lhs.derived(),rhs); } - + template<typename OtherDerived> const typename SparseProductReturnType<Derived,OtherDerived>::Type operator*(const MatrixBase<OtherDerived> &other) const; @@ -340,7 +351,7 @@ template<typename Derived> class SparseMatrixBase const SparseInnerVectorSet<Derived,1> col(int j) const; SparseInnerVectorSet<Derived,1> innerVector(int outer); const SparseInnerVectorSet<Derived,1> innerVector(int outer) const; - + // set of sub-vectors SparseInnerVectorSet<Derived,Dynamic> subrows(int start, int size); const SparseInnerVectorSet<Derived,Dynamic> subrows(int start, int size) const; @@ -432,10 +443,16 @@ template<typename Derived> class SparseMatrixBase for (int j=0; j<outerSize(); ++j) { for (typename Derived::InnerIterator i(derived(),j); i; ++i) + { if(IsRowMajor) - res.coeffRef(j,i.index()) = i.value(); + std::cerr << i.row() << "," << i.col() << " == " << j << "," << i.index() << "\n"; else - res.coeffRef(i.index(),j) = i.value(); + std::cerr << i.row() << "," << i.col() << " == " << i.index() << "," << j << "\n"; +// if(IsRowMajor) + res.coeffRef(i.row(),i.col()) = i.value(); +// else +// res.coeffRef(i.index(),j) = i.value(); + } } return res; } diff --git a/Eigen/src/Sparse/SparseRedux.h b/Eigen/src/Sparse/SparseRedux.h index f0d370548..e223c37ec 100644 --- a/Eigen/src/Sparse/SparseRedux.h +++ b/Eigen/src/Sparse/SparseRedux.h @@ -37,4 +37,20 @@ SparseMatrixBase<Derived>::sum() const return res; } +template<typename _Scalar, int _Options> +typename ei_traits<SparseMatrix<_Scalar,_Options> >::Scalar +SparseMatrix<_Scalar,_Options>::sum() const +{ + ei_assert(rows()>0 && cols()>0 && "you are using a non initialized matrix"); + return Matrix<Scalar,1,Dynamic>::Map(m_data.value(0), m_data.size()).sum(); +} + +template<typename _Scalar, int _Options> +typename ei_traits<SparseVector<_Scalar,_Options> >::Scalar +SparseVector<_Scalar,_Options>::sum() const +{ + ei_assert(rows()>0 && cols()>0 && "you are using a non initialized matrix"); + return Matrix<Scalar,1,Dynamic>::Map(m_data.value(0), m_data.size()).sum(); +} + #endif // EIGEN_SPARSEREDUX_H diff --git a/Eigen/src/Sparse/SparseTranspose.h b/Eigen/src/Sparse/SparseTranspose.h index 89a14d707..879c377a9 100644 --- a/Eigen/src/Sparse/SparseTranspose.h +++ b/Eigen/src/Sparse/SparseTranspose.h @@ -66,20 +66,26 @@ template<typename MatrixType> class SparseTranspose template<typename MatrixType> class SparseTranspose<MatrixType>::InnerIterator : public MatrixType::InnerIterator { + typedef typename MatrixType::InnerIterator Base; public: EIGEN_STRONG_INLINE InnerIterator(const SparseTranspose& trans, int outer) - : MatrixType::InnerIterator(trans.m_matrix, outer) + : Base(trans.m_matrix, outer) {} + inline int row() const { return Base::col(); } + inline int col() const { return Base::row(); } }; template<typename MatrixType> class SparseTranspose<MatrixType>::ReverseInnerIterator : public MatrixType::ReverseInnerIterator { + typedef typename MatrixType::ReverseInnerIterator Base; public: EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTranspose& xpr, int outer) - : MatrixType::ReverseInnerIterator(xpr.m_matrix, outer) + : Base(xpr.m_matrix, outer) {} + inline int row() const { return Base::col(); } + inline int col() const { return Base::row(); } }; #endif // EIGEN_SPARSETRANSPOSE_H diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h index ac815e698..05356942b 100644 --- a/Eigen/src/Sparse/SparseVector.h +++ b/Eigen/src/Sparse/SparseVector.h @@ -34,26 +34,26 @@ * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. * */ -template<typename _Scalar, int _Flags> -struct ei_traits<SparseVector<_Scalar, _Flags> > +template<typename _Scalar, int _Options> +struct ei_traits<SparseVector<_Scalar, _Options> > { typedef _Scalar Scalar; enum { - IsColVector = _Flags & RowMajorBit ? 0 : 1, + IsColVector = _Options & RowMajorBit ? 0 : 1, RowsAtCompileTime = IsColVector ? Dynamic : 1, ColsAtCompileTime = IsColVector ? 1 : Dynamic, MaxRowsAtCompileTime = RowsAtCompileTime, MaxColsAtCompileTime = ColsAtCompileTime, - Flags = SparseBit | _Flags, + Flags = SparseBit | _Options, CoeffReadCost = NumTraits<Scalar>::ReadCost, SupportedAccessPatterns = InnerRandomAccessPattern }; }; -template<typename _Scalar, int _Flags> +template<typename _Scalar, int _Options> class SparseVector - : public SparseMatrixBase<SparseVector<_Scalar, _Flags> > + : public SparseMatrixBase<SparseVector<_Scalar, _Options> > { public: EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector) @@ -124,7 +124,7 @@ class SparseVector { ei_assert(outer==0); } - + inline Scalar& insertBack(int outer, int inner) { ei_assert(outer==0); @@ -183,7 +183,7 @@ class SparseVector m_data.append(0, i); return m_data.value(m_data.size()-1); } - + /** \deprecated use insert(int,int) */ EIGEN_DEPRECATED Scalar& fillrand(int r, int c) { @@ -196,11 +196,11 @@ class SparseVector { return insert(i); } - + /** \deprecated use finalize() */ EIGEN_DEPRECATED void endFill() {} inline void finalize() {} - + void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) { m_data.prune(reference,epsilon); @@ -357,10 +357,13 @@ class SparseVector /** Destructor */ inline ~SparseVector() {} + + /** Overloaded for performance */ + Scalar sum() const; }; -template<typename Scalar, int _Flags> -class SparseVector<Scalar,_Flags>::InnerIterator +template<typename Scalar, int _Options> +class SparseVector<Scalar,_Options>::InnerIterator { public: InnerIterator(const SparseVector& vec, int outer=0) |