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