diff --git a/Eigen/src/PardisoSupport/PardisoSupport.h b/Eigen/src/PardisoSupport/PardisoSupport.h old mode 100644 new mode 100755 index 7ab2e3e6b..234e3213b --- a/Eigen/src/PardisoSupport/PardisoSupport.h +++ b/Eigen/src/PardisoSupport/PardisoSupport.h @@ -54,7 +54,7 @@ namespace internal template<> struct pardiso_run_selector { - typedef long long int IndexTypeType; + typedef long long int IndexType; static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x) { @@ -93,19 +93,19 @@ namespace internal typedef typename _MatrixType::StorageIndex StorageIndex; }; -} +} // end namespace internal template -class PardisoImpl : public SparseSolveBase +class PardisoImpl : public SparseSolverBase { protected: - typedef SparseSolveBase Base; + typedef SparseSolverBase Base; using Base::derived; using Base::m_isInitialized; typedef internal::pardiso_traits Traits; public: - using base::_solve_impl; + using Base::_solve_impl; typedef typename Traits::MatrixType MatrixType; typedef typename Traits::Scalar Scalar; @@ -173,16 +173,17 @@ class PardisoImpl : public SparseSolveBase Derived& compute(const MatrixType& matrix); - template - bool _solve_impl(const MatrixBase &b, MatrixBase& x) const; + template + void _solve_impl(const MatrixBase &b, MatrixBase &dest) const; protected: void pardisoRelease() { if(m_isInitialized) // Factorization ran at least once { - internal::pardiso_run_selector::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0, - m_iparm.data(), m_msglvl, 0, 0); + internal::pardiso_run_selector::run(m_pt, 1, 1, m_type, -1, m_size,0, 0, 0, m_perm.data(), 0, + m_iparm.data(), m_msglvl, NULL, NULL); + m_isInitialized = false; } } @@ -217,12 +218,14 @@ class PardisoImpl : public SparseSolveBase m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0; m_iparm[34] = 1; // C indexing m_iparm[59] = 1; // Automatic switch between In-Core and Out-of-Core modes + + memset(m_pt, 0, sizeof(m_pt)); } protected: // cached data to reduce reallocation, etc. - void manageErrorCode(Index error) + void manageErrorCode(Index error) const { switch(error) { @@ -239,7 +242,7 @@ class PardisoImpl : public SparseSolveBase } mutable SparseMatrixType m_matrix; - ComputationInfo m_info; + mutable ComputationInfo m_info; bool m_analysisIsOk, m_factorizationIsOk; Index m_type, m_msglvl; mutable void *m_pt[64]; @@ -256,7 +259,6 @@ Derived& PardisoImpl::compute(const MatrixType& a) eigen_assert(a.rows() == a.cols()); pardisoRelease(); - memset(m_pt, 0, sizeof(m_pt)); m_perm.setZero(m_size); derived().getMatrix(a); @@ -279,7 +281,6 @@ Derived& PardisoImpl::analyzePattern(const MatrixType& a) eigen_assert(m_size == a.cols()); pardisoRelease(); - memset(m_pt, 0, sizeof(m_pt)); m_perm.setZero(m_size); derived().getMatrix(a); @@ -313,12 +314,15 @@ Derived& PardisoImpl::factorize(const MatrixType& a) return derived(); } -template +template template -bool PardisoImpl::_solve_impl(const MatrixBase &b, MatrixBase& x) const +void PardisoImpl::_solve_impl(const MatrixBase &b, MatrixBase& x) const { if(m_iparm[0] == 0) // Factorization was not computed - return false; + { + m_info = InvalidInput; + return; + } //Index n = m_matrix.rows(); Index nrhs = Index(b.cols()); @@ -353,7 +357,7 @@ bool PardisoImpl::_solve_impl(const MatrixBase &b, MatrixBase class PardisoLU : public PardisoImpl< PardisoLU > { protected: - typedef PardisoImpl< PardisoLU > Base; + typedef PardisoImpl Base; typedef typename Base::Scalar Scalar; typedef typename Base::RealScalar RealScalar; using Base::pardisoInit; @@ -401,6 +405,7 @@ class PardisoLU : public PardisoImpl< PardisoLU > void getMatrix(const MatrixType& matrix) { m_matrix = matrix; + m_matrix.makeCompressed(); } }; @@ -424,7 +429,6 @@ class PardisoLLT : public PardisoImpl< PardisoLLT > protected: typedef PardisoImpl< PardisoLLT > Base; typedef typename Base::Scalar Scalar; - typedef typename Base::StorageIndex StorageIndex; typedef typename Base::RealScalar RealScalar; using Base::pardisoInit; using Base::m_matrix; @@ -432,9 +436,9 @@ class PardisoLLT : public PardisoImpl< PardisoLLT > public: + typedef typename Base::StorageIndex StorageIndex; enum { UpLo = _UpLo }; using Base::compute; - using Base::solve; PardisoLLT() : Base() @@ -457,6 +461,7 @@ class PardisoLLT : public PardisoImpl< PardisoLLT > PermutationMatrix p_null; m_matrix.resize(matrix.rows(), matrix.cols()); m_matrix.template selfadjointView() = matrix.template selfadjointView().twistedBy(p_null); + m_matrix.makeCompressed(); } }; @@ -482,7 +487,6 @@ class PardisoLDLT : public PardisoImpl< PardisoLDLT > protected: typedef PardisoImpl< PardisoLDLT > Base; typedef typename Base::Scalar Scalar; - typedef typename Base::StorageIndex StorageIndex; typedef typename Base::RealScalar RealScalar; using Base::pardisoInit; using Base::m_matrix; @@ -490,8 +494,8 @@ class PardisoLDLT : public PardisoImpl< PardisoLDLT > public: + typedef typename Base::StorageIndex StorageIndex; using Base::compute; - using Base::solve; enum { UpLo = Options&(Upper|Lower) }; PardisoLDLT() @@ -513,6 +517,7 @@ class PardisoLDLT : public PardisoImpl< PardisoLDLT > PermutationMatrix p_null; m_matrix.resize(matrix.rows(), matrix.cols()); m_matrix.template selfadjointView() = matrix.template selfadjointView().twistedBy(p_null); + m_matrix.makeCompressed(); } };