1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/* 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Copyright (c) 2011, Intel Corporation. All rights reserved. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Redistribution and use in source and binary forms, with or without modification, 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath are permitted provided that the following conditions are met: 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Redistributions of source code must retain the above copyright notice, this 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath list of conditions and the following disclaimer. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Redistributions in binary form must reproduce the above copyright notice, 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath this list of conditions and the following disclaimer in the documentation 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath and/or other materials provided with the distribution. 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Neither the name of Intel Corporation nor the names of its contributors may 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath be used to endorse or promote products derived from this software without 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath specific prior written permission. 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ******************************************************************************** 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Content : Eigen bindings to Intel(R) MKL PARDISO 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ******************************************************************************** 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath*/ 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_PARDISOSUPPORT_H 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_PARDISOSUPPORT_H 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> class PardisoLU; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int Options=Upper> class PardisoLLT; 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int Options=Upper> class PardisoLDLT; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename IndexType> 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct pardiso_run_selector 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, 472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x) 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang IndexType error = 0; 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error); 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return error; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<> 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct pardiso_run_selector<long long int> 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef long long int IndexType; 582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, 592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x) 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang IndexType error = 0; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error); 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return error; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<class Pardiso> struct pardiso_traits; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename _MatrixType> 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct pardiso_traits< PardisoLU<_MatrixType> > 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Scalar Scalar; 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::RealScalar RealScalar; 752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename _MatrixType::StorageIndex StorageIndex; 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename _MatrixType, int Options> 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct pardiso_traits< PardisoLLT<_MatrixType, Options> > 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Scalar Scalar; 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::RealScalar RealScalar; 842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename _MatrixType::StorageIndex StorageIndex; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename _MatrixType, int Options> 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct pardiso_traits< PardisoLDLT<_MatrixType, Options> > 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Scalar Scalar; 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::RealScalar RealScalar; 932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename _MatrixType::StorageIndex StorageIndex; 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace internal 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<class Derived> 992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass PardisoImpl : public SparseSolverBase<Derived> 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang protected: 1022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseSolverBase<Derived> Base; 1032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::derived; 1042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::m_isInitialized; 1052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::pardiso_traits<Derived> Traits; 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 1082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::_solve_impl; 1092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Traits::MatrixType MatrixType; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Traits::Scalar Scalar; 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Traits::RealScalar RealScalar; 1132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename Traits::StorageIndex StorageIndex; 1142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseMatrix<Scalar,RowMajor,StorageIndex> SparseMatrixType; 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> VectorType; 1162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<StorageIndex, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 1172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<StorageIndex, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 1182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Array<StorageIndex,64,1,DontAlign> ParameterType; 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 1202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ScalarIsComplex = NumTraits<Scalar>::IsComplex, 1212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = Dynamic, 1222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxColsAtCompileTime = Dynamic 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PardisoImpl() 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert((sizeof(StorageIndex) >= sizeof(_INTEGER_t) && sizeof(StorageIndex) <= 8) && "Non-supported index type"); 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm.setZero(); 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_msglvl = 0; // No output 1302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_isInitialized = false; 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ~PardisoImpl() 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoRelease(); 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_size; } 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_size; } 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c NumericalIssue if the matrix appears to be negative. 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_isInitialized && "Decomposition is not initialized."); 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \warning for advanced usage only. 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a reference to the parameter array controlling PARDISO. 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * See the PARDISO manual to know how to use it. */ 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ParameterType& pardisoParameterArray() 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_iparm; 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Performs a symbolic decomposition on the sparcity of \a matrix. 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function is particularly useful when solving for several problems having the same structure. 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Derived& analyzePattern(const MatrixType& matrix); 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Performs a numeric decomposition of \a matrix 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Derived& factorize(const MatrixType& matrix); 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Derived& compute(const MatrixType& matrix); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Rhs,typename Dest> 1792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void pardisoRelease() 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_isInitialized) // Factorization ran at least once 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, -1, internal::convert_index<StorageIndex>(m_size),0, 0, 0, m_perm.data(), 0, 1872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_iparm.data(), m_msglvl, NULL, NULL); 1882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_isInitialized = false; 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void pardisoInit(int type) 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_type = type; 195615d816d068b4d0f5e8df601930b5f160bf7eda1Tim Murray bool symmetric = std::abs(m_type) < 10; 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[0] = 1; // No solver default 1972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_iparm[1] = 2; // use Metis for the ordering 1982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_iparm[2] = 0; // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??) 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[3] = 0; // No iterative-direct algorithm 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[4] = 0; // No user fill-in reducing permutation 2012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_iparm[5] = 0; // Write solution into x, b is left unchanged 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[6] = 0; // Not in use 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[7] = 2; // Max numbers of iterative refinement steps 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[8] = 0; // Not in use 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[9] = 13; // Perturb the pivot elements with 1E-13 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[11] = 0; // Not in use 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric). 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Try m_iparm[12] = 1 in case of inappropriate accuracy 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[13] = 0; // Output: Number of perturbed pivots 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[14] = 0; // Not in use 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[15] = 0; // Not in use 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[16] = 0; // Not in use 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[18] = -1; // Output: Mflops for LU factorization 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[19] = 0; // Output: Numbers of CG Iterations 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[20] = 0; // 1x1 pivoting 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[26] = 0; // No matrix checker 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0; 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[34] = 1; // C indexing 2222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_iparm[36] = 0; // CSR 2232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_iparm[59] = 0; // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core 2242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang memset(m_pt, 0, sizeof(m_pt)); 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // cached data to reduce reallocation, etc. 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void manageErrorCode(Index error) const 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath switch(error) 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case 0: 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case -4: 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case -7: 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NumericalIssue; 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath default: 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = InvalidInput; 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable SparseMatrixType m_matrix; 2482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang mutable ComputationInfo m_info; 2492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang bool m_analysisIsOk, m_factorizationIsOk; 2502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex m_type, m_msglvl; 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable void *m_pt[64]; 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable ParameterType m_iparm; 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable IntColVectorType m_perm; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index m_size; 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<class Derived> 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathDerived& PardisoImpl<Derived>::compute(const MatrixType& a) 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_size = a.rows(); 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(a.rows() == a.cols()); 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoRelease(); 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_perm.setZero(m_size); 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derived().getMatrix(a); 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index error; 2692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size), 2702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 2712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath manageErrorCode(error); 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = true; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = true; 2752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_isInitialized = true; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<class Derived> 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathDerived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a) 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_size = a.rows(); 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_size == a.cols()); 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoRelease(); 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_perm.setZero(m_size); 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derived().getMatrix(a); 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index error; 2902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, internal::convert_index<StorageIndex>(m_size), 2912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 2922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath manageErrorCode(error); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = true; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = false; 2972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_isInitialized = true; 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<class Derived> 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathDerived& PardisoImpl<Derived>::factorize(const MatrixType& a) 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_size == a.rows() && m_size == a.cols()); 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derived().getMatrix(a); 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index error; 3102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, internal::convert_index<StorageIndex>(m_size), 3112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 3122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath manageErrorCode(error); 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = true; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<class Derived> 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename BDerived,typename XDerived> 3212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid PardisoImpl<Derived>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_iparm[0] == 0) // Factorization was not computed 3242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_info = InvalidInput; 3262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return; 3272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //Index n = m_matrix.rows(); 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nrhs = Index(b.cols()); 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_size==b.rows()); 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported"); 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported"); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(((nrhs == 1) || b.outerStride() == b.rows())); 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// switch (transposed) { 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// case SvNoTrans : m_iparm[11] = 0 ; break; 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// case SvTranspose : m_iparm[11] = 2 ; break; 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// case SvAdjoint : m_iparm[11] = 1 ; break; 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// default: 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n"; 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// m_iparm[11] = 0; 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// } 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data()); 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp; 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Pardiso cannot solve in-place 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rhs_ptr == x.derived().data()) 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp = b; 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs_ptr = tmp.data(); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index error; 3572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, internal::convert_index<StorageIndex>(m_size), 3582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 3592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl, 3602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rhs_ptr, x.derived().data()); 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang manageErrorCode(error); 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PardisoSupport_Module 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PardisoLU 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A sparse direct LU factorization and solver based on the PARDISO library 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * using the Intel MKL PARDISO library. The sparse matrix A must be squared and invertible. 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse. 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 3742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set: 3752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \code solver.pardisoParameterArray()[59] = 1; \endcode 3762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 3792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \implsparsesolverconcept 3802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 3812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa \ref TutorialSparseSolverConcept, class SparseLU 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PardisoLU : public PardisoImpl< PardisoLU<MatrixType> > 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 3872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef PardisoImpl<PardisoLU> Base; 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::Scalar Scalar; 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::RealScalar RealScalar; 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::pardisoInit; 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_matrix; 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath friend class PardisoImpl< PardisoLU<MatrixType> >; 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::compute; 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::solve; 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PardisoLU() 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base() 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoInit(Base::ScalarIsComplex ? 13 : 11); 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit PardisoLU(const MatrixType& matrix) 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base() 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoInit(Base::ScalarIsComplex ? 13 : 11); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void getMatrix(const MatrixType& matrix) 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix = matrix; 4152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.makeCompressed(); 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PardisoSupport_Module 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PardisoLLT 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class allows to solve for A.X = B sparse linear problems via a LL^T Cholesky factorization 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * using the Intel MKL PARDISO library. The sparse matrix A must be selfajoint and positive definite. 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse. 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 4272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set: 4282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \code solver.pardisoParameterArray()[59] = 1; \endcode 4292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam UpLo can be any bitwise combination of Upper, Lower. The default is Upper, meaning only the upper triangular part has to be used. 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Upper|Lower can be used to tell both triangular parts can be used as input. 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 4342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \implsparsesolverconcept 4352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 4362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa \ref TutorialSparseSolverConcept, class SimplicialLLT 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> > 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base; 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::Scalar Scalar; 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::RealScalar RealScalar; 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::pardisoInit; 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_matrix; 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >; 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename Base::StorageIndex StorageIndex; 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::compute; 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PardisoLLT() 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base() 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoInit(Base::ScalarIsComplex ? 4 : 2); 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit PardisoLLT(const MatrixType& matrix) 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base() 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoInit(Base::ScalarIsComplex ? 4 : 2); 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void getMatrix(const MatrixType& matrix) 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // PARDISO supports only upper, row-major matrices 4732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null; 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.resize(matrix.rows(), matrix.cols()); 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); 4762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.makeCompressed(); 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PardisoSupport_Module 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PardisoLDLT 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class allows to solve for A.X = B sparse linear problems via a LDL^T Cholesky factorization 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * using the Intel MKL PARDISO library. The sparse matrix A is assumed to be selfajoint and positive definite. 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * For complex matrices, A can also be symmetric only, see the \a Options template parameter. 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse. 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 4892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set: 4902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \code solver.pardisoParameterArray()[59] = 1; \endcode 4912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam Options can be any bitwise combination of Upper, Lower, and Symmetric. The default is Upper, meaning only the upper triangular part has to be used. 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Symmetric can be used for symmetric, non-selfadjoint complex matrices, the default being to assume a selfadjoint matrix. 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Upper|Lower can be used to tell both triangular parts can be used as input. 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 4972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \implsparsesolverconcept 4982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 4992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa \ref TutorialSparseSolverConcept, class SimplicialLDLT 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int Options> 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base; 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::Scalar Scalar; 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::RealScalar RealScalar; 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::pardisoInit; 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_matrix; 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >; 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename Base::StorageIndex StorageIndex; 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::compute; 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = Options&(Upper|Lower) }; 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PardisoLDLT() 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base() 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit PardisoLDLT(const MatrixType& matrix) 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base() 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void getMatrix(const MatrixType& matrix) 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // PARDISO supports only upper, row-major matrices 5342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null; 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.resize(matrix.rows(), matrix.cols()); 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); 5372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.makeCompressed(); 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_PARDISOSUPPORT_H 544