diff options
author | Gael Guennebaud <g.gael@free.fr> | 2012-02-27 13:23:21 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2012-02-27 13:23:21 +0100 |
commit | 122f28626c86a2e5965a58c6a69a31f7ada4e149 (patch) | |
tree | f2d98689988c92c243d31667861ce558ae1fab9f | |
parent | b240a3fad9d0147aeb287a1ee74238558081a2ec (diff) |
fix and clean Pardiso solver and s/PARDISOSupport/PardisoSupport
-rw-r--r-- | Eigen/PardisoSupport (renamed from Eigen/PARDISOSupport) | 6 | ||||
-rw-r--r-- | Eigen/src/PardisoSupport/CMakeLists.txt (renamed from Eigen/src/PARDISOSupport/CMakeLists.txt) | 0 | ||||
-rw-r--r-- | Eigen/src/PardisoSupport/PardisoSupport.h (renamed from Eigen/src/PARDISOSupport/PARDISOSupport.h) | 294 | ||||
-rw-r--r-- | test/pardiso_support.cpp | 14 |
4 files changed, 201 insertions, 113 deletions
diff --git a/Eigen/PARDISOSupport b/Eigen/PardisoSupport index 3d079b18b..d11dac171 100644 --- a/Eigen/PARDISOSupport +++ b/Eigen/PardisoSupport @@ -12,16 +12,16 @@ namespace Eigen { /** \ingroup Support_modules - * \defgroup PARDISOSupport_Module PARDISOSupport module + * \defgroup PardisoSupport_Module PardisoSupport module * * This module brings support for the Intel(R) MKL PARDISO direct sparse solvers * * \code - * #include <Eigen/PARDISOSupport> + * #include <Eigen/PardisoSupport> * \endcode */ -#include "src/PARDISOSupport/PARDISOSupport.h" +#include "src/PardisoSupport/PardisoSupport.h" } // namespace Eigen diff --git a/Eigen/src/PARDISOSupport/CMakeLists.txt b/Eigen/src/PardisoSupport/CMakeLists.txt index a097ab401..a097ab401 100644 --- a/Eigen/src/PARDISOSupport/CMakeLists.txt +++ b/Eigen/src/PardisoSupport/CMakeLists.txt diff --git a/Eigen/src/PARDISOSupport/PARDISOSupport.h b/Eigen/src/PardisoSupport/PardisoSupport.h index ca3b66ca4..3904e6418 100644 --- a/Eigen/src/PARDISOSupport/PARDISOSupport.h +++ b/Eigen/src/PardisoSupport/PardisoSupport.h @@ -73,8 +73,8 @@ namespace internal typedef typename _MatrixType::Index Index; }; - template<typename _MatrixType> - struct pardiso_traits< PardisoLLT<_MatrixType> > + template<typename _MatrixType, int Options> + struct pardiso_traits< PardisoLLT<_MatrixType, Options> > { typedef _MatrixType MatrixType; typedef typename _MatrixType::Scalar Scalar; @@ -82,13 +82,13 @@ namespace internal typedef typename _MatrixType::Index Index; }; - template<typename _MatrixType> - struct pardiso_traits< PardisoLDLT<_MatrixType> > + template<typename _MatrixType, int Options> + struct pardiso_traits< PardisoLDLT<_MatrixType, Options> > { typedef _MatrixType MatrixType; typedef typename _MatrixType::Scalar Scalar; typedef typename _MatrixType::RealScalar RealScalar; - typedef typename _MatrixType::Index Index; + typedef typename _MatrixType::Index Index; }; } @@ -96,11 +96,13 @@ namespace internal template<class Derived> class PardisoImpl { + typedef internal::pardiso_traits<Derived> Traits; public: - typedef typename internal::pardiso_traits<Derived>::MatrixType MatrixType; - typedef typename internal::pardiso_traits<Derived>::Scalar Scalar; - typedef typename internal::pardiso_traits<Derived>::RealScalar RealScalar; - typedef typename internal::pardiso_traits<Derived>::Index Index; + typedef typename Traits::MatrixType MatrixType; + typedef typename Traits::Scalar Scalar; + typedef typename Traits::RealScalar RealScalar; + typedef typename Traits::Index Index; + typedef SparseMatrix<Scalar,RowMajor,Index> SparseMatrixType; typedef Matrix<Scalar,Dynamic,1> VectorType; typedef Matrix<Index, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; typedef Matrix<Index, MatrixType::RowsAtCompileTime, 1> IntColVectorType; @@ -112,7 +114,7 @@ class PardisoImpl { eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type"); m_iparm.setZero(); - m_msglvl = 0; /* No output */ + m_msglvl = 0; // No output m_initialized = false; } @@ -121,8 +123,8 @@ class PardisoImpl pardisoRelease(); } - inline Index cols() const { return m_matrix.cols(); } - inline Index rows() const { return m_matrix.rows(); } + inline Index cols() const { return m_size; } + inline Index rows() const { return m_size; } /** \brief Reports whether previous computation was successful. * @@ -142,8 +144,25 @@ class PardisoImpl { return m_param; } + + /** Performs a symbolic decomposition on the sparcity of \a matrix. + * + * This function is particularly useful when solving for several problems having the same structure. + * + * \sa factorize() + */ + Derived& analyzePattern(const MatrixType& matrix); + + /** Performs a numeric decomposition of \a matrix + * + * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. + * + * \sa analyzePattern() + */ + Derived& factorize(const MatrixType& matrix); Derived& compute(const MatrixType& matrix); + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -188,19 +207,22 @@ class PardisoImpl template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const { - eigen_assert(m_matrix.rows()==b.rows()); + 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(); - Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols); + // 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.leftCols(actualCols) = b.middleCols(k,actualCols); - tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols)); - dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView(); + 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(); } } @@ -209,9 +231,8 @@ class PardisoImpl { if(m_initialized) // Factorization ran at least once { - internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_matrix.rows(), NULL, NULL, NULL, m_perm.data(), 0, - m_iparm.data(), m_msglvl, NULL, NULL); - m_iparm.setZero(); + internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0, + m_iparm.data(), m_msglvl, 0, 0); } } @@ -219,106 +240,142 @@ class PardisoImpl { m_type = type; bool symmetric = abs(m_type) < 10; - m_iparm[0] = 1; /* No solver default */ - m_iparm[1] = 3; // use Metis for the ordering - /* Numbers of processors, value of OMP_NUM_THREADS */ - m_iparm[2] = 1; - m_iparm[3] = 0; /* No iterative-direct algorithm */ - m_iparm[4] = 0; /* No user fill-in reducing permutation */ - m_iparm[5] = 0; /* Write solution into x */ - m_iparm[6] = 0; /* Not in use */ - m_iparm[7] = 2; /* Max numbers of iterative refinement steps */ - m_iparm[8] = 0; /* Not in use */ - m_iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */ - m_iparm[10] = symmetric ? 0 : 1; /* Use nonsymmetric permutation and scaling MPS */ - m_iparm[11] = 0; /* Not in use */ - m_iparm[12] = symmetric ? 0 : 1; /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try m_iparm[12] = 1 in case of inappropriate accuracy */ - m_iparm[13] = 0; /* Output: Number of perturbed pivots */ - m_iparm[14] = 0; /* Not in use */ - m_iparm[15] = 0; /* Not in use */ - m_iparm[16] = 0; /* Not in use */ - m_iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */ - m_iparm[18] = -1; /* Output: Mflops for LU factorization */ - m_iparm[19] = 0; /* Output: Numbers of CG Iterations */ - m_iparm[20] = 0; /* 1x1 pivoting */ - m_iparm[26] = 0; /* No matrix checker */ + m_iparm[0] = 1; // No solver default + m_iparm[1] = 3; // use Metis for the ordering + m_iparm[2] = 1; // Numbers of processors, value of OMP_NUM_THREADS + m_iparm[3] = 0; // No iterative-direct algorithm + m_iparm[4] = 0; // No user fill-in reducing permutation + m_iparm[5] = 0; // Write solution into x + m_iparm[6] = 0; // Not in use + m_iparm[7] = 2; // Max numbers of iterative refinement steps + m_iparm[8] = 0; // Not in use + m_iparm[9] = 13; // Perturb the pivot elements with 1E-13 + m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS + m_iparm[11] = 0; // Not in use + m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric). + // Try m_iparm[12] = 1 in case of inappropriate accuracy + m_iparm[13] = 0; // Output: Number of perturbed pivots + m_iparm[14] = 0; // Not in use + m_iparm[15] = 0; // Not in use + m_iparm[16] = 0; // Not in use + m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU + m_iparm[18] = -1; // Output: Mflops for LU factorization + m_iparm[19] = 0; // Output: Numbers of CG Iterations + + m_iparm[20] = 0; // 1x1 pivoting + m_iparm[26] = 0; // No matrix checker m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0; - m_iparm[34] = 0; /* Fortran indexing */ - m_iparm[59] = 1; /* Automatic switch between In-Core and Out-of-Core modes */ + m_iparm[34] = 1; // C indexing + m_iparm[59] = 1; // Automatic switch between In-Core and Out-of-Core modes } protected: // cached data to reduce reallocation, etc. + + void manageErrorCode(Index error) + { + switch(error) + { + case 0: + m_info = Success; + break; + case -4: + case -7: + m_info = NumericalIssue; + break; + default: + m_info = InvalidInput; + } + } + mutable SparseMatrixType m_matrix; ComputationInfo m_info; - bool m_initialized, m_succeeded; + bool m_initialized, m_analysisIsOk, m_factorizationIsOk; Index m_type, m_msglvl; mutable void *m_pt[64]; mutable Array<Index,64,1> m_iparm; - mutable SparseMatrix<Scalar, RowMajor> m_matrix; mutable IntColVectorType m_perm; + Index m_size; }; template<class Derived> Derived& PardisoImpl<Derived>::compute(const MatrixType& a) { - Index n = a.rows(), i; + m_size = a.rows(); eigen_assert(a.rows() == a.cols()); pardisoRelease(); memset(m_pt, 0, sizeof(m_pt)); + m_perm.setZero(m_size); + derived().getMatrix(a); + + Index error; + error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size, + m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), + m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); + + manageErrorCode(error); + m_analysisIsOk = true; + m_factorizationIsOk = true; m_initialized = true; + return derived(); +} - bool symmetric = abs(m_type) < 10; - m_iparm[10] = symmetric ? 0 : 1; /* Use nonsymmetric permutation and scaling MPS */ - m_iparm[12] = symmetric ? 0 : 1; /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try m_iparm[12] = 1 in case of inappropriate accuracy */ - - m_perm.resize(n); - m_matrix = a; +template<class Derived> +Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a) +{ + m_size = a.rows(); + eigen_assert(m_size == a.cols()); - /* Convert to Fortran-style indexing */ - for(i = 0; i <= m_matrix.rows(); ++i) - ++m_matrix.outerIndexPtr()[i]; - for(i = 0; i < m_matrix.nonZeros(); ++i) - ++m_matrix.innerIndexPtr()[i]; + pardisoRelease(); + memset(m_pt, 0, sizeof(m_pt)); + m_perm.setZero(m_size); + derived().getMatrix(a); + + Index error; + error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size, + m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), + m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); + + manageErrorCode(error); + m_analysisIsOk = true; + m_factorizationIsOk = false; + m_initialized = true; + return derived(); +} - Index error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, n, - m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), - m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); +template<class Derived> +Derived& PardisoImpl<Derived>::factorize(const MatrixType& a) +{ + eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); + eigen_assert(m_size == a.rows() && m_size == a.cols()); + + derived().getMatrix(a); - switch(error) - { - case 0: - m_succeeded = true; - m_info = Success; - return derived(); - case -4: - case -7: - m_info = NumericalIssue; - break; - default: - m_info = InvalidInput; - } - m_succeeded = false; + Index error; + error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size, + m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), + m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); + + manageErrorCode(error); + m_factorizationIsOk = true; return derived(); } template<class Base> template<typename BDerived,typename XDerived> -bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, - MatrixBase<XDerived>& x) const +bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const { if(m_iparm[0] == 0) // Factorization was not computed return false; - Index n = m_matrix.rows(); + //Index n = m_matrix.rows(); Index nrhs = b.cols(); - eigen_assert(n==b.rows()); + eigen_assert(m_size==b.rows()); eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported"); eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported"); eigen_assert(((nrhs == 1) || b.outerStride() == b.rows())); - //x.derived().resizeLike(b); // switch (transposed) { // case SvNoTrans : m_iparm[11] = 0 ; break; @@ -329,15 +386,27 @@ bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, // m_iparm[11] = 0; // } - Index error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, n, - m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), - m_perm.data(), nrhs, m_iparm.data(), m_msglvl, const_cast<Scalar*>(&b(0, 0)), &x(0, 0)); + Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data()); + Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp; + + // Pardiso cannot solve in-place + if(rhs_ptr == x.derived().data()) + { + tmp = b; + rhs_ptr = tmp.data(); + } + + Index error; + error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size, + m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), + m_perm.data(), nrhs, m_iparm.data(), m_msglvl, + rhs_ptr, x.derived().data()); return error==0; } -/** \ingroup PARDISOSupport_Module +/** \ingroup PardisoSupport_Module * \class PardisoLU * \brief A sparse direct LU factorization and solver based on the PARDISO library * @@ -357,6 +426,8 @@ class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> > typedef typename Base::Scalar Scalar; typedef typename Base::RealScalar RealScalar; using Base::pardisoInit; + using Base::m_matrix; + friend class PardisoImpl< PardisoLU<MatrixType> >; public: @@ -375,9 +446,14 @@ class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> > pardisoInit(Base::ScalarIsComplex ? 13 : 11); compute(matrix); } + protected: + void getMatrix(const MatrixType& matrix) + { + m_matrix = matrix; + } }; -/** \ingroup PARDISOSupport_Module +/** \ingroup PardisoSupport_Module * \class PardisoLLT * \brief A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library * @@ -395,10 +471,13 @@ template<typename MatrixType, int _UpLo> class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> > { protected: - typedef PardisoImpl< PardisoLLT<MatrixType> > Base; + typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base; typedef typename Base::Scalar Scalar; + typedef typename Base::Index Index; typedef typename Base::RealScalar RealScalar; using Base::pardisoInit; + using Base::m_matrix; + friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >; public: @@ -418,11 +497,21 @@ class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> > pardisoInit(Base::ScalarIsComplex ? 4 : 2); compute(matrix); } + + protected: + + void getMatrix(const MatrixType& matrix) + { + // PARDISO supports only upper, row-major matrices + PermutationMatrix<Dynamic,Dynamic,Index> p_null; + m_matrix.resize(matrix.rows(), matrix.cols()); + m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); + } }; -/** \ingroup PARDISOSupport_Module +/** \ingroup PardisoSupport_Module * \class PardisoLDLT - * \brief A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library + * \brief A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library * * This class allows to solve for A.X = B sparse linear problems via a LDL^T Cholesky factorization * using the Intel MKL PARDISO library. The sparse matrix A is assumed to be selfajoint and positive definite. @@ -440,11 +529,13 @@ template<typename MatrixType, int Options> class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > { protected: - typedef PardisoImpl< PardisoLDLT<MatrixType> > Base; + typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base; typedef typename Base::Scalar Scalar; typedef typename Base::Index Index; typedef typename Base::RealScalar RealScalar; using Base::pardisoInit; + using Base::m_matrix; + friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >; public: @@ -459,26 +550,19 @@ class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > } PardisoLDLT(const MatrixType& matrix) - : Base(flags) + : Base() { pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); - compute(matrix, hermitian); + compute(matrix); } - - void compute(const MatrixType& matrix) + + void getMatrix(const MatrixType& matrix) { - if(Options&Upper==0) - { - // PARDISO supports only upper, row-major matrices - PermutationMatrix<Dynamic,Dynamic,Index> P(0); - SparseMatrix<Scalar,RowMajor> tmp(matrix.rows(), matrix.cols()); - tmp.template selfadjointView<Upper>() = matrix.template selfadjointView<Lower>().twistedBy(P); - Base::compute(tmp); - } - else - Base::compute(matrix); + // PARDISO supports only upper, row-major matrices + PermutationMatrix<Dynamic,Dynamic,Index> p_null; + m_matrix.resize(matrix.rows(), matrix.cols()); + m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); } - }; namespace internal { diff --git a/test/pardiso_support.cpp b/test/pardiso_support.cpp index 316c608d0..11cb98e10 100644 --- a/test/pardiso_support.cpp +++ b/test/pardiso_support.cpp @@ -3,16 +3,20 @@ */ #include "sparse_solver.h" -#include <Eigen/PARDISOSupport> +#include <Eigen/PardisoSupport> template<typename T> void test_pardiso_T() { - PardisoLLT < SparseMatrix<T, RowMajor> > pardiso_llt; - PardisoLDLT< SparseMatrix<T, RowMajor> > pardiso_ldlt; + PardisoLLT < SparseMatrix<T, RowMajor>, Lower> pardiso_llt_lower; + PardisoLLT < SparseMatrix<T, RowMajor>, Upper> pardiso_llt_upper; + PardisoLDLT < SparseMatrix<T, RowMajor>, Lower> pardiso_ldlt_lower; + PardisoLDLT < SparseMatrix<T, RowMajor>, Upper> pardiso_ldlt_upper; PardisoLU < SparseMatrix<T, RowMajor> > pardiso_lu; - check_sparse_spd_solving(pardiso_llt); - check_sparse_spd_solving(pardiso_ldlt); + check_sparse_spd_solving(pardiso_llt_lower); + check_sparse_spd_solving(pardiso_llt_upper); + check_sparse_spd_solving(pardiso_ldlt_lower); + check_sparse_spd_solving(pardiso_ldlt_upper); check_sparse_square_solving(pardiso_lu); } |