sparse module: * add row(i), col(i) functions * add prune() function to remove small coefficients
diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h index 350630c..6005568 100644 --- a/Eigen/src/Core/util/StaticAssert.h +++ b/Eigen/src/Core/util/StaticAssert.h
@@ -72,6 +72,7 @@ INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION, YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES, + THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES, INVALID_MATRIX_TEMPLATE_PARAMETERS }; };
diff --git a/Eigen/src/Sparse/CompressedStorage.h b/Eigen/src/Sparse/CompressedStorage.h index d66cb5c..4e8074b 100644 --- a/Eigen/src/Sparse/CompressedStorage.h +++ b/Eigen/src/Sparse/CompressedStorage.h
@@ -31,6 +31,7 @@ template<typename Scalar> class CompressedStorage { + typedef typename NumTraits<Scalar>::Real RealScalar; public: CompressedStorage() : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) @@ -183,6 +184,22 @@ } return m_values[id]; } + + void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) + { + int k = 0; + int n = size(); + for (int i=0; i<n; ++i) + { + if (!ei_isMuchSmallerThan(value(i), reference, epsilon)) + { + value(k) = value(i); + index(k) = index(i); + ++k; + } + } + resize(k,0); + } protected:
diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h index 30e54a9..5e7ce93 100644 --- a/Eigen/src/Sparse/DynamicSparseMatrix.h +++ b/Eigen/src/Sparse/DynamicSparseMatrix.h
@@ -74,7 +74,7 @@ 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; } @@ -102,8 +102,6 @@ return m_data[outer].atWithInsertion(inner); } - public: - class InnerIterator; inline void setZero() @@ -176,6 +174,12 @@ /** Does nothing. Provided for compatibility with SparseMatrix. */ inline void endFill() {} + + void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) + { + for (int j=0; j<outerSize(); ++j) + m_data[j].prune(reference,epsilon); + } /** Resize the matrix without preserving the data (the matrix is set to zero) */
diff --git a/Eigen/src/Sparse/SparseBlock.h b/Eigen/src/Sparse/SparseBlock.h index 48be375..cbec36e 100644 --- a/Eigen/src/Sparse/SparseBlock.h +++ b/Eigen/src/Sparse/SparseBlock.h
@@ -78,6 +78,40 @@ } }; +/** \returns the i-th row of the matrix \c *this. For row-major matrix only. */ +template<typename Derived> +SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i) +{ + EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES); + return innerVector(i); +} + +/** \returns the i-th row of the matrix \c *this. For row-major matrix only. + * (read-only version) */ +template<typename Derived> +const SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i) const +{ + EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES); + return innerVector(i); +} + +/** \returns the i-th column of the matrix \c *this. For column-major matrix only. */ +template<typename Derived> +SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i) +{ + EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES); + return innerVector(i); +} + +/** \returns the i-th column of the matrix \c *this. For column-major matrix only. + * (read-only version) */ +template<typename Derived> +const SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i) const +{ + EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES); + return innerVector(i); +} + /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this * is col-major (resp. row-major). */
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 2ce4664..e0c56be 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -69,11 +69,11 @@ int* m_outerIndex; CompressedStorage<Scalar> m_data; - public: inline int rows() const { return IsRowMajor ? m_outerSize : m_innerSize; } inline int cols() const { return IsRowMajor ? m_innerSize : m_outerSize; } + inline int innerSize() const { return m_innerSize; } inline int outerSize() const { return m_outerSize; } inline int innerNonZeros(int j) const { return m_outerIndex[j+1]-m_outerIndex[j]; } @@ -230,6 +230,28 @@ ++i; } } + + void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) + { + int k = 0; + for (int j=0; j<m_outerSize; ++j) + { + int previousStart = m_outerIndex[j]; + m_outerIndex[j] = k; + int end = m_outerIndex[j+1]; + for (int i=previousStart; i<end; ++i) + { + if (!ei_isMuchSmallerThan(m_data.value(i), reference, epsilon)) + { + m_data.value(k) = m_data.value(i); + m_data.index(k) = m_data.index(i); + ++k; + } + } + } + m_outerIndex[m_outerSize] = k; + m_data.resize(k,0); + } void resize(int rows, int cols) {
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 32011af..46c0c54 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -327,6 +327,10 @@ // void transposeInPlace(); const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; } + SparseInnerVector<Derived> row(int i); + const SparseInnerVector<Derived> row(int i) const; + SparseInnerVector<Derived> col(int j); + const SparseInnerVector<Derived> col(int j) const; SparseInnerVector<Derived> innerVector(int outer); const SparseInnerVector<Derived> innerVector(int outer) const;
diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h index ff9d27a..9af5653 100644 --- a/Eigen/src/Sparse/SparseVector.h +++ b/Eigen/src/Sparse/SparseVector.h
@@ -67,20 +67,19 @@ CompressedStorage<Scalar> m_data; int m_size; - public: - inline int rows() const { return IsColVector ? m_size : 1; } - inline int cols() const { return IsColVector ? 1 : m_size; } - inline int innerSize() const { return m_size; } - inline int outerSize() const { return 1; } - inline int innerNonZeros(int j) const { ei_assert(j==0); return m_size; } + EIGEN_STRONG_INLINE int rows() const { return IsColVector ? m_size : 1; } + EIGEN_STRONG_INLINE int cols() const { return IsColVector ? 1 : m_size; } + EIGEN_STRONG_INLINE int innerSize() const { return m_size; } + EIGEN_STRONG_INLINE int outerSize() const { return 1; } + EIGEN_STRONG_INLINE int innerNonZeros(int j) const { ei_assert(j==0); return m_size; } - inline const Scalar* _valuePtr() const { return &m_data.value(0); } - inline Scalar* _valuePtr() { return &m_data.value(0); } + EIGEN_STRONG_INLINE const Scalar* _valuePtr() const { return &m_data.value(0); } + EIGEN_STRONG_INLINE Scalar* _valuePtr() { return &m_data.value(0); } - inline const int* _innerIndexPtr() const { return &m_data.index(0); } - inline int* _innerIndexPtr() { return &m_data.index(0); } + EIGEN_STRONG_INLINE const int* _innerIndexPtr() const { return &m_data.index(0); } + EIGEN_STRONG_INLINE int* _innerIndexPtr() { return &m_data.index(0); } inline Scalar coeff(int row, int col) const { @@ -145,6 +144,11 @@ m_data.value(id+1) = 0; return m_data.value(id+1); } + + void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) + { + m_data.prune(reference,epsilon); + } void resize(int newSize) {
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index addd40f..5ae038e 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp
@@ -323,14 +323,49 @@ VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b); } + + // test prune + { + SparseMatrixType m2(rows, rows); + DenseMatrix refM2(rows, rows); + refM2.setZero(); + int countFalseNonZero = 0; + int countTrueNonZero = 0; + m2.startFill(); + for (int j=0; j<m2.outerSize(); ++j) + for (int i=0; i<m2.innerSize(); ++i) + { + float x = ei_random<float>(0,1); + if (x<0.1) + { + // do nothing + } + else if (x<0.5) + { + countFalseNonZero++; + m2.fill(i,j) = Scalar(0); + } + else + { + countTrueNonZero++; + m2.fill(i,j) = refM2(i,j) = Scalar(1); + } + } + m2.endFill(); + VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros()); + VERIFY_IS_APPROX(m2, refM2); + m2.prune(1); + VERIFY(countTrueNonZero==m2.nonZeros()); + VERIFY_IS_APPROX(m2, refM2); + } } void test_sparse_basic() { for(int i = 0; i < g_repeat; i++) { -// CALL_SUBTEST( sparse_basic(SparseMatrix<double>(8, 8)) ); -// CALL_SUBTEST( sparse_basic(SparseMatrix<std::complex<double> >(16, 16)) ); -// CALL_SUBTEST( sparse_basic(SparseMatrix<double>(33, 33)) ); + CALL_SUBTEST( sparse_basic(SparseMatrix<double>(8, 8)) ); + CALL_SUBTEST( sparse_basic(SparseMatrix<std::complex<double> >(16, 16)) ); + CALL_SUBTEST( sparse_basic(SparseMatrix<double>(33, 33)) ); CALL_SUBTEST( sparse_basic(DynamicSparseMatrix<double>(8, 8)) ); }