diff options
-rw-r--r-- | Eigen/src/PaStiXSupport/PaStiXSupport.h | 23 | ||||
-rw-r--r-- | Eigen/src/PardisoSupport/PardisoSupport.h | 25 | ||||
-rw-r--r-- | Eigen/src/SparseCholesky/SimplicialCholesky.h | 23 | ||||
-rw-r--r-- | Eigen/src/SuperLUSupport/SuperLUSupport.h | 18 | ||||
-rw-r--r-- | Eigen/src/UmfPackSupport/UmfPackSupport.h | 21 | ||||
-rw-r--r-- | Eigen/src/misc/SparseSolve.h | 17 | ||||
-rw-r--r-- | cmake/FindSuperLU.cmake | 3 | ||||
-rw-r--r-- | test/sparse_solver.h | 11 |
8 files changed, 50 insertions, 91 deletions
diff --git a/Eigen/src/PaStiXSupport/PaStiXSupport.h b/Eigen/src/PaStiXSupport/PaStiXSupport.h index 82e137c64..a955287d1 100644 --- a/Eigen/src/PaStiXSupport/PaStiXSupport.h +++ b/Eigen/src/PaStiXSupport/PaStiXSupport.h @@ -157,27 +157,6 @@ class PastixBase : internal::noncopyable template<typename Rhs,typename Dest> bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const; - /** \internal */ - template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> - void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &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()"); - eigen_assert(rows()==b.rows()); - - // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. - static const int NbColsAtOnce = 1; - int rhsCols = b.cols(); - int size = b.rows(); - Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols); - for(int k=0; k<rhsCols; k+=NbColsAtOnce) - { - int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); - tmp.leftCols(actualCols) = b.middleCols(k,actualCols); - tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols)); - dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView(); - } - } - Derived& derived() { return *static_cast<Derived*>(this); @@ -731,7 +710,7 @@ struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve_sparse(rhs(),dst); + this->defaultEvalTo(dst); } }; diff --git a/Eigen/src/PardisoSupport/PardisoSupport.h b/Eigen/src/PardisoSupport/PardisoSupport.h index d623bf518..1c48f0df7 100644 --- a/Eigen/src/PardisoSupport/PardisoSupport.h +++ b/Eigen/src/PardisoSupport/PardisoSupport.h @@ -206,29 +206,6 @@ class PardisoImpl template<typename BDerived, typename XDerived> bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const; - /** \internal */ - template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> - void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const - { - eigen_assert(m_size==b.rows()); - - // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. - static const int NbColsAtOnce = 4; - int rhsCols = b.cols(); - int size = b.rows(); - // Pardiso cannot solve in-place, - // so we need two temporaries - Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_rhs(size,rhsCols); - Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_res(size,rhsCols); - for(int k=0; k<rhsCols; k+=NbColsAtOnce) - { - int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); - tmp_rhs.leftCols(actualCols) = b.middleCols(k,actualCols); - tmp_res.leftCols(actualCols) = derived().solve(tmp_rhs.leftCols(actualCols)); - dest.middleCols(k,actualCols) = tmp_res.leftCols(actualCols).sparseView(); - } - } - protected: void pardisoRelease() { @@ -604,7 +581,7 @@ struct sparse_solve_retval<PardisoImpl<Derived>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec().derived()._solve_sparse(rhs(),dst); + this->defaultEvalTo(dst); } }; diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index 51f6fe9ef..59bddb1e4 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -215,27 +215,6 @@ class SimplicialCholeskyBase : internal::noncopyable dest = m_Pinv * dest; } - /** \internal */ - template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> - void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &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()"); - eigen_assert(m_matrix.rows()==b.rows()); - - // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. - static const int NbColsAtOnce = 4; - int rhsCols = b.cols(); - int size = b.rows(); - Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols); - for(int k=0; k<rhsCols; k+=NbColsAtOnce) - { - int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); - tmp.leftCols(actualCols) = b.middleCols(k,actualCols); - tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols)); - dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView(); - } - } - #endif // EIGEN_PARSED_BY_DOXYGEN protected: @@ -864,7 +843,7 @@ struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec().derived()._solve_sparse(rhs(),dst); + this->defaultEvalTo(dst); } }; diff --git a/Eigen/src/SuperLUSupport/SuperLUSupport.h b/Eigen/src/SuperLUSupport/SuperLUSupport.h index cd6c4b91f..3034c7af5 100644 --- a/Eigen/src/SuperLUSupport/SuperLUSupport.h +++ b/Eigen/src/SuperLUSupport/SuperLUSupport.h @@ -353,14 +353,14 @@ class SuperLUBase : internal::noncopyable * * \sa compute() */ -// template<typename Rhs> -// inline const internal::sparse_solve_retval<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const -// { -// eigen_assert(m_isInitialized && "SuperLU is not initialized."); -// eigen_assert(rows()==b.rows() -// && "SuperLU::solve(): invalid number of rows of the right hand side matrix b"); -// return internal::sparse_solve_retval<SuperLU, Rhs>(*this, b.derived()); -// } + template<typename Rhs> + inline const internal::sparse_solve_retval<SuperLUBase, Rhs> solve(const SparseMatrixBase<Rhs>& b) const + { + eigen_assert(m_isInitialized && "SuperLU is not initialized."); + eigen_assert(rows()==b.rows() + && "SuperLU::solve(): invalid number of rows of the right hand side matrix b"); + return internal::sparse_solve_retval<SuperLUBase, Rhs>(*this, b.derived()); + } /** Performs a symbolic decomposition on the sparcity of \a matrix. * @@ -1015,7 +1015,7 @@ struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec().derived()._solve(rhs(),dst); + this->defaultEvalTo(dst); } }; diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h index 22d049089..d85b8be85 100644 --- a/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -215,14 +215,14 @@ class UmfPackLU : internal::noncopyable * * \sa compute() */ -// template<typename Rhs> -// inline const internal::sparse_solve_retval<UmfPAckLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const -// { -// eigen_assert(m_isInitialized && "UmfPAckLU is not initialized."); -// eigen_assert(rows()==b.rows() -// && "UmfPAckLU::solve(): invalid number of rows of the right hand side matrix b"); -// return internal::sparse_solve_retval<UmfPAckLU, Rhs>(*this, b.derived()); -// } + template<typename Rhs> + inline const internal::sparse_solve_retval<UmfPackLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const + { + eigen_assert(m_isInitialized && "UmfPackLU is not initialized."); + eigen_assert(rows()==b.rows() + && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b"); + return internal::sparse_solve_retval<UmfPackLU, Rhs>(*this, b.derived()); + } /** Performs a symbolic decomposition on the sparcity of \a matrix. * @@ -381,7 +381,8 @@ bool UmfPackLU<MatrixType>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDe const int rhsCols = b.cols(); eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet"); eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet"); - + eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve"); + int errorCode; for (int j=0; j<rhsCols; ++j) { @@ -420,7 +421,7 @@ struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + this->defaultEvalTo(dst); } }; diff --git a/Eigen/src/misc/SparseSolve.h b/Eigen/src/misc/SparseSolve.h index 272c4a479..244bb8ec7 100644 --- a/Eigen/src/misc/SparseSolve.h +++ b/Eigen/src/misc/SparseSolve.h @@ -47,6 +47,23 @@ template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_b } protected: + template<typename DestScalar, int DestOptions, typename DestIndex> + inline void defaultEvalTo(SparseMatrix<DestScalar,DestOptions,DestIndex>& dst) const + { + // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. + static const int NbColsAtOnce = 4; + int rhsCols = m_rhs.cols(); + int size = m_rhs.rows(); + Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols); + Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,rhsCols); + for(int k=0; k<rhsCols; k+=NbColsAtOnce) + { + int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); + tmp.leftCols(actualCols) = m_rhs.middleCols(k,actualCols); + tmpX.leftCols(actualCols) = m_dec.solve(tmp.leftCols(actualCols)); + dst.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView(); + } + } const DecompositionType& m_dec; typename Rhs::Nested m_rhs; }; diff --git a/cmake/FindSuperLU.cmake b/cmake/FindSuperLU.cmake index ca72b4498..8a3df3666 100644 --- a/cmake/FindSuperLU.cmake +++ b/cmake/FindSuperLU.cmake @@ -14,9 +14,10 @@ find_path(SUPERLU_INCLUDES ${INCLUDE_INSTALL_DIR} PATH_SUFFIXES superlu + SRC ) -find_library(SUPERLU_LIBRARIES superlu PATHS $ENV{SUPERLUDIR} ${LIB_INSTALL_DIR}) +find_library(SUPERLU_LIBRARIES superlu PATHS $ENV{SUPERLUDIR} ${LIB_INSTALL_DIR} PATH_SUFFIXES lib) include(FindPackageHandleStandardArgs) find_package_handle_standard_args(SUPERLU DEFAULT_MSG diff --git a/test/sparse_solver.h b/test/sparse_solver.h index 73d92874c..5a1be67e7 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -37,7 +37,6 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); VERIFY(x.isApprox(refX,test_precision<Scalar>())); - x.setZero(); // test the analyze/factorize API solver.analyzePattern(A); @@ -258,6 +257,7 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver) { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; + typedef SparseMatrix<Scalar,ColMajor> SpMat; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; @@ -267,12 +267,17 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver) DenseMatrix dA; int size = generate_sparse_square_problem(solver, A, dA); - DenseVector b = DenseVector::Random(size); - DenseMatrix dB = DenseMatrix::Random(size,rhsCols); A.makeCompressed(); + DenseVector b = DenseVector::Random(size); + DenseMatrix dB(size,rhsCols); + SpMat B(size,rhsCols); + double density = (std::max)(8./(size*rhsCols), 0.1); + initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); + B.makeCompressed(); for (int i = 0; i < g_repeat; i++) { check_sparse_solving(solver, A, b, dA, b); check_sparse_solving(solver, A, dB, dA, dB); + check_sparse_solving(solver, A, B, dA, dB); } // First, get the folder |