From 241e5ee3e7a306feb7b70f7f57b94c7e045abb67 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 27 Oct 2010 14:31:23 +0200 Subject: [PATCH] add the possibility to solve for sparse rhs with Cholmod --- Eigen/src/Sparse/DynamicSparseMatrix.h | 19 ++-- Eigen/src/Sparse/SparseMatrix.h | 9 ++ Eigen/src/Sparse/SparseMatrixBase.h | 7 ++ Eigen/src/Sparse/SparseVector.h | 4 + .../Eigen/src/SparseExtra/CholmodSupport.h | 106 +++++++++++++++++- unsupported/test/sparse_llt.cpp | 43 +++++-- 6 files changed, 171 insertions(+), 17 deletions(-) diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h index 2c85a1d96..88ad85f84 100644 --- a/Eigen/src/Sparse/DynamicSparseMatrix.h +++ b/Eigen/src/Sparse/DynamicSparseMatrix.h @@ -44,8 +44,8 @@ */ namespace internal { -template -struct traits > +template +struct traits > { typedef _Scalar Scalar; typedef _Index Index; @@ -56,16 +56,16 @@ struct traits > ColsAtCompileTime = Dynamic, MaxRowsAtCompileTime = Dynamic, MaxColsAtCompileTime = Dynamic, - Flags = _Flags | NestByRefBit | LvalueBit, + Flags = _Options | NestByRefBit | LvalueBit, CoeffReadCost = NumTraits::ReadCost, SupportedAccessPatterns = OuterRandomAccessPattern }; }; } -template +template class DynamicSparseMatrix - : public SparseMatrixBase > + : public SparseMatrixBase > { public: EIGEN_SPARSE_PUBLIC_INTERFACE(DynamicSparseMatrix) @@ -74,6 +74,9 @@ class DynamicSparseMatrix // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=) typedef MappedSparseMatrix Map; using Base::IsRowMajor; + enum { + Options = _Options + }; protected: @@ -325,10 +328,10 @@ class DynamicSparseMatrix EIGEN_DEPRECATED void endFill() {} }; -template -class DynamicSparseMatrix::InnerIterator : public SparseVector::InnerIterator +template +class DynamicSparseMatrix::InnerIterator : public SparseVector::InnerIterator { - typedef typename SparseVector::InnerIterator Base; + typedef typename SparseVector::InnerIterator Base; public: InnerIterator(const DynamicSparseMatrix& mat, Index outer) : Base(mat.m_data[outer]), m_outer(outer) diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index c8b75a96c..9a6e0b763 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -78,6 +78,9 @@ class SparseMatrix typedef MappedSparseMatrix Map; using Base::IsRowMajor; typedef CompressedStorage Storage; + enum { + Options = _Options + }; protected: @@ -436,6 +439,12 @@ class SparseMatrix { return Base::operator=(product); } + + template + EIGEN_STRONG_INLINE SparseMatrix& operator=(const ReturnByValue& func) + { + return Base::operator=(func); + } #endif template diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 8750de12c..ba6f64ab7 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -183,6 +183,13 @@ template class SparseMatrixBase : public EigenBase this->operator=(other); return derived(); } + + template + Derived& operator=(const ReturnByValue& other) + { + other.evalTo(derived()); + return derived(); + } template diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h index be1ad4c50..c640eadff 100644 --- a/Eigen/src/Sparse/SparseVector.h +++ b/Eigen/src/Sparse/SparseVector.h @@ -72,6 +72,10 @@ class SparseVector typedef SparseMatrixBase SparseBase; enum { IsColVector = internal::traits::IsColVector }; + + enum { + Options = _Options + }; CompressedStorage m_data; Index m_size; diff --git a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h index 08f01a6e6..b0f1c8e42 100644 --- a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h +++ b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h @@ -26,7 +26,57 @@ #define EIGEN_CHOLMODSUPPORT_H namespace internal { + +template struct sparse_solve_retval_base; +template struct sparse_solve_retval; +template +struct traits > +{ + typedef typename DecompositionType::MatrixType MatrixType; + typedef SparseMatrix ReturnType; +}; + +template struct sparse_solve_retval_base + : public ReturnByValue > +{ + typedef typename remove_all::type RhsNestedCleaned; + typedef _DecompositionType DecompositionType; + typedef ReturnByValue Base; + typedef typename Base::Index Index; + + sparse_solve_retval_base(const DecompositionType& dec, const Rhs& rhs) + : m_dec(dec), m_rhs(rhs) + {} + + inline Index rows() const { return m_dec.cols(); } + inline Index cols() const { return m_rhs.cols(); } + inline const DecompositionType& dec() const { return m_dec; } + inline const RhsNestedCleaned& rhs() const { return m_rhs; } + + template inline void evalTo(Dest& dst) const + { + static_cast*>(this)->evalTo(dst); + } + + protected: + const DecompositionType& m_dec; + const typename Rhs::Nested m_rhs; +}; + +#define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType,Rhs) \ + typedef typename DecompositionType::MatrixType MatrixType; \ + typedef typename MatrixType::Scalar Scalar; \ + typedef typename MatrixType::RealScalar RealScalar; \ + typedef typename MatrixType::Index Index; \ + typedef Eigen::internal::sparse_solve_retval_base Base; \ + using Base::dec; \ + using Base::rhs; \ + using Base::rows; \ + using Base::cols; \ + sparse_solve_retval(const DecompositionType& dec, const Rhs& rhs) \ + : Base(dec, rhs) {} + template void cholmod_configure_matrix(CholmodType& mat) { @@ -94,6 +144,13 @@ cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat) return res; } +template +const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat) +{ + cholmod_sparse res = viewAsCholmod(mat.const_cast_derived()); + return res; +} + /** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix. * The data are not copied but shared. */ template @@ -247,6 +304,20 @@ class CholmodDecomposition return internal::solve_retval(*this, b.derived()); } + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. + * + * \sa compute() + */ + template + inline const internal::sparse_solve_retval + solve(const SparseMatrixBase& b) const + { + eigen_assert(m_isInitialized && "LLT is not initialized."); + eigen_assert(rows()==b.rows() + && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b"); + return internal::sparse_solve_retval(*this, b.derived()); + } + /** Performs a symbolic decomposition on the sparcity of \a matrix. * * This function is particularly useful when solving for several problems having the same structure. @@ -305,10 +376,30 @@ class CholmodDecomposition { this->m_info = NumericalIssue; } - dest = Matrix::Map(reinterpret_cast(x_cd->x),b.rows()); + // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.) + dest = Matrix::Map(reinterpret_cast(x_cd->x),b.rows(),b.cols()); cholmod_free_dense(&x_cd, &m_cholmod); } + /** \internal */ + template + void _solve(const SparseMatrix &b, SparseMatrix &dest) const + { + eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); + const Index size = m_cholmodFactor->n; + eigen_assert(size==b.rows()); + + // note: cs stands for Cholmod Sparse + cholmod_sparse b_cs = viewAsCholmod(b); + cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod); + if(!x_cs) + { + this->m_info = NumericalIssue; + } + // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.) + dest = viewAsEigen(*x_cs); + cholmod_free_sparse(&x_cs, &m_cholmod); + } #endif // EIGEN_PARSED_BY_DOXYGEN protected: @@ -335,6 +426,19 @@ struct solve_retval, Rhs> } }; +template +struct sparse_solve_retval, Rhs> + : sparse_solve_retval_base, Rhs> +{ + typedef CholmodDecomposition<_MatrixType,_UpLo> Dec; + EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) + + template void evalTo(Dest& dst) const + { + dec()._solve(rhs(),dst); + } +}; + } #endif // EIGEN_CHOLMODSUPPORT_H diff --git a/unsupported/test/sparse_llt.cpp b/unsupported/test/sparse_llt.cpp index 21bd36d35..7b1cad4c5 100644 --- a/unsupported/test/sparse_llt.cpp +++ b/unsupported/test/sparse_llt.cpp @@ -40,19 +40,21 @@ template void sparse_llt(int rows, int cols) DenseMatrix refMat2(rows, cols); DenseVector b = DenseVector::Random(cols); - DenseVector refX(cols), x(cols); + DenseVector ref_x(cols), x(cols); + DenseMatrix B = DenseMatrix::Random(rows,cols); + DenseMatrix ref_X(rows,cols), X(rows,cols); initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0); for(int i=0; i().llt().solve(b); + ref_x = refMat2.template selfadjointView().llt().solve(b); if (!NumTraits::IsComplex) { x = b; SparseLLT > (m2).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision()) && "LLT: default"); + VERIFY(ref_x.isApprox(x,test_precision()) && "LLT: default"); } #ifdef EIGEN_CHOLMOD_SUPPORT @@ -62,14 +64,14 @@ template void sparse_llt(int rows, int cols) SparseMatrix m3 = m2.adjoint()*m2; DenseMatrix refMat3 = refMat2.adjoint()*refMat2; - refX = refMat3.template selfadjointView().llt().solve(b); + ref_x = refMat3.template selfadjointView().llt().solve(b); x = b; SparseLLT, Cholmod>(m3).solveInPlace(x); VERIFY((m3*x).isApprox(b,test_precision()) && "LLT legacy: cholmod solveInPlace"); x = SparseLLT, Cholmod>(m3).solve(b); - VERIFY(refX.isApprox(x,test_precision()) && "LLT legacy: cholmod solve"); + VERIFY(ref_x.isApprox(x,test_precision()) && "LLT legacy: cholmod solve"); } // new API @@ -78,13 +80,38 @@ template void sparse_llt(int rows, int cols) SparseMatrix m3 = m2.adjoint()*m2; DenseMatrix refMat3 = refMat2.adjoint()*refMat2; - refX = refMat3.template selfadjointView().llt().solve(b); + // with a single vector as the rhs + ref_x = refMat3.template selfadjointView().llt().solve(b); x = CholmodDecomposition, Lower>(m3).solve(b); - VERIFY(refX.isApprox(x,test_precision()) && "LLT: cholmod solve"); + VERIFY(ref_x.isApprox(x,test_precision()) && "LLT: cholmod solve, single dense rhs"); x = CholmodDecomposition, Upper>(m3).solve(b); - VERIFY(refX.isApprox(x,test_precision()) && "LLT: cholmod solve"); + VERIFY(ref_x.isApprox(x,test_precision()) && "LLT: cholmod solve, single dense rhs"); + + + // with multiple rhs + ref_X = refMat3.template selfadjointView().llt().solve(B); + + X = CholmodDecomposition, Lower>(m3).solve(B); + VERIFY(ref_X.isApprox(X,test_precision()) && "LLT: cholmod solve, multiple dense rhs"); + + X = CholmodDecomposition, Upper>(m3).solve(B); + VERIFY(ref_X.isApprox(X,test_precision()) && "LLT: cholmod solve, multiple dense rhs"); + + + // with a sparse rhs + SparseMatrix spB(rows,cols), spX(rows,cols); + B.diagonal().array() += 1; + spB = B.sparseView(0.5,1); + + ref_X = refMat3.template selfadjointView().llt().solve(DenseMatrix(spB)); + + spX = CholmodDecomposition, Lower>(m3).solve(spB); + VERIFY(ref_X.isApprox(spX.toDense(),test_precision()) && "LLT: cholmod solve, multiple sparse rhs"); + + spX = CholmodDecomposition, Upper>(m3).solve(spB); + VERIFY(ref_X.isApprox(spX.toDense(),test_precision()) && "LLT: cholmod solve, multiple sparse rhs"); } #endif