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 { runpardiso_run_selector46 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 { 127 eigen_assert((sizeof(StorageIndex) >= sizeof(_INTEGER_t) && sizeof(StorageIndex) <= 8) && "Non-supported index type"); 128 m_iparm.setZero(); 129 m_msglvl = 0; // No output 130 m_isInitialized = false; 131 } 132 133 ~PardisoImpl() 134 { 135 pardisoRelease(); 136 } 137 138 inline Index cols() const { return m_size; } 139 inline Index rows() const { return m_size; } 140 141 /** \brief Reports whether previous computation was successful. 142 * 143 * \returns \c Success if computation was succesful, 144 * \c NumericalIssue if the matrix appears to be negative. 145 */ 146 ComputationInfo info() const 147 { 148 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 149 return m_info; 150 } 151 152 /** \warning for advanced usage only. 153 * \returns a reference to the parameter array controlling PARDISO. 154 * See the PARDISO manual to know how to use it. */ 155 ParameterType& pardisoParameterArray() 156 { 157 return m_iparm; 158 } 159 160 /** Performs a symbolic decomposition on the sparcity of \a matrix. 161 * 162 * This function is particularly useful when solving for several problems having the same structure. 163 * 164 * \sa factorize() 165 */ 166 Derived& analyzePattern(const MatrixType& matrix); 167 168 /** Performs a numeric decomposition of \a matrix 169 * 170 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 171 * 172 * \sa analyzePattern() 173 */ 174 Derived& factorize(const MatrixType& matrix); 175 176 Derived& compute(const MatrixType& matrix); 177 178 template<typename Rhs,typename Dest> 179 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; 180 181 protected: 182 void pardisoRelease() 183 { 184 if(m_isInitialized) // Factorization ran at least once 185 { 186 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, 187 m_iparm.data(), m_msglvl, NULL, NULL); 188 m_isInitialized = false; 189 } 190 } 191 192 void pardisoInit(int type) 193 { 194 m_type = type; 195 EIGEN_USING_STD_MATH(abs); 196 bool symmetric = 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 typedef typename Base::Scalar Scalar; 390 typedef typename Base::RealScalar RealScalar; 391 using Base::pardisoInit; 392 using Base::m_matrix; 393 friend class PardisoImpl< PardisoLU<MatrixType> >; 394 395 public: 396 397 using Base::compute; 398 using Base::solve; 399 400 PardisoLU() 401 : Base() 402 { 403 pardisoInit(Base::ScalarIsComplex ? 13 : 11); 404 } 405 406 explicit PardisoLU(const MatrixType& matrix) 407 : Base() 408 { 409 pardisoInit(Base::ScalarIsComplex ? 13 : 11); 410 compute(matrix); 411 } 412 protected: 413 void getMatrix(const MatrixType& matrix) 414 { 415 m_matrix = matrix; 416 m_matrix.makeCompressed(); 417 } 418 }; 419 420 /** \ingroup PardisoSupport_Module 421 * \class PardisoLLT 422 * \brief A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library 423 * 424 * This class allows to solve for A.X = B sparse linear problems via a LL^T Cholesky factorization 425 * using the Intel MKL PARDISO library. The sparse matrix A must be selfajoint and positive definite. 426 * The vectors or matrices X and B can be either dense or sparse. 427 * 428 * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set: 429 * \code solver.pardisoParameterArray()[59] = 1; \endcode 430 * 431 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 432 * \tparam UpLo can be any bitwise combination of Upper, Lower. The default is Upper, meaning only the upper triangular part has to be used. 433 * Upper|Lower can be used to tell both triangular parts can be used as input. 434 * 435 * \implsparsesolverconcept 436 * 437 * \sa \ref TutorialSparseSolverConcept, class SimplicialLLT 438 */ 439 template<typename MatrixType, int _UpLo> 440 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> > 441 { 442 protected: 443 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base; 444 typedef typename Base::Scalar Scalar; 445 typedef typename Base::RealScalar RealScalar; 446 using Base::pardisoInit; 447 using Base::m_matrix; 448 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >; 449 450 public: 451 452 typedef typename Base::StorageIndex StorageIndex; 453 enum { UpLo = _UpLo }; 454 using Base::compute; 455 456 PardisoLLT() 457 : Base() 458 { 459 pardisoInit(Base::ScalarIsComplex ? 4 : 2); 460 } 461 462 explicit PardisoLLT(const MatrixType& matrix) 463 : Base() 464 { 465 pardisoInit(Base::ScalarIsComplex ? 4 : 2); 466 compute(matrix); 467 } 468 469 protected: 470 471 void getMatrix(const MatrixType& matrix) 472 { 473 // PARDISO supports only upper, row-major matrices 474 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null; 475 m_matrix.resize(matrix.rows(), matrix.cols()); 476 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); 477 m_matrix.makeCompressed(); 478 } 479 }; 480 481 /** \ingroup PardisoSupport_Module 482 * \class PardisoLDLT 483 * \brief A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library 484 * 485 * This class allows to solve for A.X = B sparse linear problems via a LDL^T Cholesky factorization 486 * using the Intel MKL PARDISO library. The sparse matrix A is assumed to be selfajoint and positive definite. 487 * For complex matrices, A can also be symmetric only, see the \a Options template parameter. 488 * The vectors or matrices X and B can be either dense or sparse. 489 * 490 * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set: 491 * \code solver.pardisoParameterArray()[59] = 1; \endcode 492 * 493 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 494 * \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. 495 * Symmetric can be used for symmetric, non-selfadjoint complex matrices, the default being to assume a selfadjoint matrix. 496 * Upper|Lower can be used to tell both triangular parts can be used as input. 497 * 498 * \implsparsesolverconcept 499 * 500 * \sa \ref TutorialSparseSolverConcept, class SimplicialLDLT 501 */ 502 template<typename MatrixType, int Options> 503 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > 504 { 505 protected: 506 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base; 507 typedef typename Base::Scalar Scalar; 508 typedef typename Base::RealScalar RealScalar; 509 using Base::pardisoInit; 510 using Base::m_matrix; 511 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >; 512 513 public: 514 515 typedef typename Base::StorageIndex StorageIndex; 516 using Base::compute; 517 enum { UpLo = Options&(Upper|Lower) }; 518 519 PardisoLDLT() 520 : Base() 521 { 522 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 523 } 524 525 explicit PardisoLDLT(const MatrixType& matrix) 526 : Base() 527 { 528 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 529 compute(matrix); 530 } 531 532 void getMatrix(const MatrixType& matrix) 533 { 534 // PARDISO supports only upper, row-major matrices 535 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null; 536 m_matrix.resize(matrix.rows(), matrix.cols()); 537 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); 538 m_matrix.makeCompressed(); 539 } 540 }; 541 542 } // end namespace Eigen 543 544 #endif // EIGEN_PARDISOSUPPORT_H 545