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 #ifndef EIGEN_USE_MKL 36 extern "C" 37 { 38 void pardiso( void *, int *, int *, int *, int *, int * , void *, int *, 39 int * , int *, int *, int *, int *, void *, void *, int *); 40 41 void pardiso_chkmatrix (int *, int *, void *, int *, int *, int *); 42 void pardiso_chkvec (int *, int *, void *, int *); 43 void pardiso_printstats(int *, int *, void *, int *, int *, int *,void *, int *); 44 45 } // extern "C" 46 #endif 47 48 namespace Eigen { 49 50 template<typename _MatrixType> class PardisoLU; 51 template<typename _MatrixType, int Options=Upper> class PardisoLLT; 52 template<typename _MatrixType, int Options=Upper> class PardisoLDLT; 53 54 namespace internal 55 { 56 template<typename IndexType> 57 struct pardiso_run_selector 58 { runpardiso_run_selector59 static IndexType run( void * pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, 60 IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x) 61 { 62 IndexType error = 0; 63 ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error); 64 return error; 65 } 66 67 68 #ifndef EIGEN_USE_MKL printstatspardiso_run_selector69 static void printstats( void * pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, 70 IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x) 71 { 72 IndexType error = 0; 73 ::pardiso_printstats(&type, &n, a, ia, ja, &nrhs, b, &error); 74 } 75 #endif 76 }; 77 78 79 #ifdef EIGEN_USE_MKL 80 template<> 81 struct pardiso_run_selector<long long int> 82 { 83 typedef long long int IndexType; // for mkl 84 static IndexType run( void * pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, 85 IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x) 86 { 87 IndexType error = 0; 88 ::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error); 89 return error; 90 } 91 }; 92 #endif 93 94 template<class Pardiso> struct pardiso_traits; 95 96 template<typename _MatrixType> 97 struct pardiso_traits< PardisoLU<_MatrixType> > 98 { 99 typedef _MatrixType MatrixType; 100 typedef typename _MatrixType::Scalar Scalar; 101 typedef typename _MatrixType::RealScalar RealScalar; 102 typedef typename _MatrixType::StorageIndex StorageIndex; 103 }; 104 105 template<typename _MatrixType, int Options> 106 struct pardiso_traits< PardisoLLT<_MatrixType, Options> > 107 { 108 typedef _MatrixType MatrixType; 109 typedef typename _MatrixType::Scalar Scalar; 110 typedef typename _MatrixType::RealScalar RealScalar; 111 typedef typename _MatrixType::StorageIndex StorageIndex; 112 }; 113 114 template<typename _MatrixType, int Options> 115 struct pardiso_traits< PardisoLDLT<_MatrixType, Options> > 116 { 117 typedef _MatrixType MatrixType; 118 typedef typename _MatrixType::Scalar Scalar; 119 typedef typename _MatrixType::RealScalar RealScalar; 120 typedef typename _MatrixType::StorageIndex StorageIndex; 121 }; 122 123 } // end namespace internal 124 125 template<class Derived> 126 class PardisoImpl : public SparseSolverBase<Derived> 127 { 128 protected: 129 typedef SparseSolverBase<Derived> Base; 130 using Base::derived; 131 using Base::m_isInitialized; 132 133 typedef internal::pardiso_traits<Derived> Traits; 134 public: 135 using Base::_solve_impl; 136 137 typedef typename Traits::MatrixType MatrixType; 138 typedef typename Traits::Scalar Scalar; 139 typedef typename Traits::RealScalar RealScalar; 140 typedef typename Traits::StorageIndex StorageIndex; 141 typedef SparseMatrix<Scalar,RowMajor,StorageIndex> SparseMatrixType; 142 typedef Matrix<Scalar,Dynamic,1> VectorType; 143 typedef Matrix<StorageIndex, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 144 typedef Matrix<StorageIndex, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 145 typedef Array<StorageIndex,64,1,DontAlign> ParameterType; 146 enum { 147 ScalarIsComplex = NumTraits<Scalar>::IsComplex, 148 ColsAtCompileTime = Dynamic, 149 MaxColsAtCompileTime = Dynamic 150 }; 151 152 PardisoImpl() 153 { 154 //eigen_assert((sizeof(StorageIndex) >= sizeof(_INTEGER_t) && sizeof(StorageIndex) <= 8) && "Non-supported index type"); 155 m_iparm.setZero(); 156 m_msglvl = 0; // No output 157 m_isInitialized = false; 158 } 159 160 ~PardisoImpl() 161 { 162 pardisoRelease(); 163 } 164 165 inline Index cols() const { return m_size; } 166 inline Index rows() const { return m_size; } 167 168 /** \brief Reports whether previous computation was successful. 169 * 170 * \returns \c Success if computation was succesful, 171 * \c NumericalIssue if the matrix appears to be negative. 172 */ 173 ComputationInfo info() const 174 { 175 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 176 return m_info; 177 } 178 179 /** \warning for advanced usage only. 180 * \returns a reference to the parameter array controlling PARDISO. 181 * See the PARDISO manual to know how to use it. */ 182 ParameterType& pardisoParameterArray() 183 { 184 return m_iparm; 185 } 186 187 // G+Smo: set indivisual parameter 188 void setParam(const int i, const int value) { m_iparm[i] = value; } 189 190 /** Performs a symbolic decomposition on the sparcity of \a matrix. 191 * 192 * This function is particularly useful when solving for several problems having the same structure. 193 * 194 * \sa factorize() 195 */ 196 Derived& analyzePattern(const MatrixType& matrix); 197 198 /** Performs a numeric decomposition of \a matrix 199 * 200 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 201 * 202 * \sa analyzePattern() 203 */ 204 Derived& factorize(const MatrixType& matrix); 205 206 Derived& compute(const MatrixType& matrix); 207 208 template<typename Rhs,typename Dest> 209 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; 210 211 protected: 212 void pardisoRelease() 213 { 214 if(m_isInitialized) // Factorization ran at least once 215 { 216 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, 217 m_iparm.data(), m_msglvl, NULL, NULL); 218 m_isInitialized = false; 219 } 220 } 221 222 void pardisoInit(int type) 223 { 224 m_type = type; 225 bool symmetric = std::abs(m_type) < 10; 226 m_iparm[0] = 1; // 0: default values 227 m_iparm[1] = // 2: use Metis for the ordering, 3: OpenMP enabled 228 #ifdef GISMO_WITH_OPENMP 229 3; 230 #else 231 2; 232 #endif 233 m_iparm[2] = 0; // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??) 234 m_iparm[3] = 0; // No iterative-direct algorithm 235 m_iparm[4] = 0; // No user fill-in reducing permutation 236 m_iparm[5] = 0; // Write solution into x, b is left unchanged 237 m_iparm[6] = 0; // Not in use (user permutation) 238 m_iparm[7] = 2; // Max numbers of iterative refinement steps 239 m_iparm[8] = 0; // Not in use 240 m_iparm[9] = 13; // Perturb the pivot elements with 1E-13 241 m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS 242 m_iparm[11] = 0; // Not in use (solve transposed matrix) 243 m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric). 244 // Try m_iparm[12] = 1 in case of inappropriate accuracy 245 m_iparm[13] = 0; // Output: Number of perturbed pivots 246 m_iparm[14] = 0; // Not in use 247 m_iparm[15] = 0; // Not in use 248 m_iparm[16] = 0; // Not in use 249 m_iparm[17] = 0; // Output: Number of nonzeros in the factor LU 250 m_iparm[18] = 0; // Output: Mflops for LU factorization 251 m_iparm[19] = 0; // Output: Numbers of CG Iterations 252 253 m_iparm[20] = 0; // 1x1 pivoting 254 m_iparm[26] = 0; // No matrix checker 255 m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0; 256 m_iparm[34] = 0; // 1: C-style indexing (MKL only) 257 m_iparm[36] = 0; // use CSR format 258 m_iparm[59] = 0; // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core 259 260 memset(m_pt, 0, sizeof(m_pt)); 261 } 262 263 protected: 264 // cached data to reduce reallocation, etc. 265 266 void manageErrorCode(Index error) const 267 { 268 switch(error) 269 { 270 case 0: 271 m_info = Success; 272 break; 273 case -4: 274 case -7: 275 m_info = NumericalIssue; 276 break; 277 default: 278 m_info = InvalidInput; 279 } 280 } 281 282 mutable SparseMatrixType m_matrix; 283 mutable ComputationInfo m_info; 284 bool m_analysisIsOk, m_factorizationIsOk; 285 StorageIndex m_type, m_msglvl; 286 mutable void *m_pt[64]; 287 mutable ParameterType m_iparm; 288 mutable IntColVectorType m_perm; 289 Index m_size; 290 291 }; 292 293 template<class Derived> 294 Derived& PardisoImpl<Derived>::compute(const MatrixType& a) 295 { 296 m_size = a.rows(); 297 eigen_assert(a.rows() == a.cols()); 298 299 pardisoRelease(); 300 m_perm.setZero(m_size); 301 derived().getMatrix(a); 302 303 Index error; 304 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size), 305 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 306 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 307 manageErrorCode(error); 308 m_analysisIsOk = true; 309 m_factorizationIsOk = true; 310 m_isInitialized = true; 311 return derived(); 312 } 313 314 template<class Derived> 315 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a) 316 { 317 m_size = a.rows(); 318 eigen_assert(m_size == a.cols()); 319 320 pardisoRelease(); 321 m_perm.setZero(m_size); 322 derived().getMatrix(a); 323 324 Index error; 325 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, internal::convert_index<StorageIndex>(m_size), 326 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 327 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 328 329 manageErrorCode(error); 330 m_analysisIsOk = true; 331 m_factorizationIsOk = false; 332 m_isInitialized = true; 333 return derived(); 334 } 335 336 template<class Derived> 337 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a) 338 { 339 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 340 eigen_assert(m_size == a.rows() && m_size == a.cols()); 341 342 derived().getMatrix(a); 343 344 Index error; 345 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, internal::convert_index<StorageIndex>(m_size), 346 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 347 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 348 349 manageErrorCode(error); 350 m_factorizationIsOk = true; 351 return derived(); 352 } 353 354 template<class Derived> 355 template<typename BDerived,typename XDerived> 356 void PardisoImpl<Derived>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const 357 { 358 if(m_iparm[0] == 0) // Factorization was not computed 359 { 360 m_info = InvalidInput; 361 return; 362 } 363 364 //Index n = m_matrix.rows(); 365 Index nrhs = Index(b.cols()); 366 eigen_assert(m_size==b.rows()); 367 eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported"); 368 eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported"); 369 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows())); 370 371 372 // switch (transposed) { 373 // case SvNoTrans : m_iparm[11] = 0 ; break; 374 // case SvTranspose : m_iparm[11] = 2 ; break; 375 // case SvAdjoint : m_iparm[11] = 1 ; break; 376 // default: 377 // //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n"; 378 // m_iparm[11] = 0; 379 // } 380 381 Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data()); 382 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp; 383 384 // Pardiso cannot solve in-place 385 if(rhs_ptr == x.derived().data()) 386 { 387 tmp = b; 388 rhs_ptr = tmp.data(); 389 } 390 391 Index error; 392 393 // Following lines check and write out details on the CSR matrix 394 /* 395 #ifndef EIGEN_USE_MKL 396 internal::pardiso_run_selector<StorageIndex>::printstats(m_pt, 1, 1, m_type, 33, internal::convert_index<StorageIndex>(m_size), 397 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 398 m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl, 399 rhs_ptr, x.derived().data()); 400 #endif 401 */ 402 403 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, internal::convert_index<StorageIndex>(m_size), 404 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 405 m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl, 406 rhs_ptr, x.derived().data()); 407 408 manageErrorCode(error); 409 } 410 411 412 /** \ingroup PardisoSupport_Module 413 * \class PardisoLU 414 * \brief A sparse direct LU factorization and solver based on the PARDISO library 415 * 416 * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization 417 * using the Intel MKL PARDISO library. The sparse matrix A must be squared and invertible. 418 * The vectors or matrices X and B can be either dense or sparse. 419 * 420 * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set: 421 * \code solver.pardisoParameterArray()[59] = 1; \endcode 422 * 423 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 424 * 425 * \implsparsesolverconcept 426 * 427 * \sa \ref TutorialSparseSolverConcept, class SparseLU 428 */ 429 template<typename MatrixType> 430 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> > 431 { 432 protected: 433 typedef PardisoImpl<PardisoLU> Base; 434 using Base::pardisoInit; 435 using Base::m_matrix; 436 friend class PardisoImpl< PardisoLU<MatrixType> >; 437 438 public: 439 440 typedef typename Base::Scalar Scalar; 441 typedef typename Base::RealScalar RealScalar; 442 443 using Base::compute; 444 using Base::solve; 445 446 PardisoLU() 447 : Base() 448 { 449 pardisoInit(Base::ScalarIsComplex ? 13 : 11); 450 } 451 452 explicit PardisoLU(const MatrixType& matrix) 453 : Base() 454 { 455 pardisoInit(Base::ScalarIsComplex ? 13 : 11); 456 compute(matrix); 457 } 458 protected: 459 void getMatrix(const MatrixType& matrix) 460 { 461 m_matrix = matrix; 462 m_matrix.makeCompressed(); 463 //1-index based conversion 464 typedef typename MatrixType::StorageIndex StorageIndex; 465 StorageIndex * s = m_matrix.innerIndexPtr(); 466 std::transform(s, s + m_matrix.nonZeros(), s, 467 internal::bind2nd_op<std::plus<StorageIndex> >(1) 468 //std::bind2nd(std::plus<StorageIndex>(), 1) 469 ); 470 s = m_matrix.outerIndexPtr(); 471 std::transform(s, s + m_matrix.outerSize()+1, s, 472 internal::bind2nd_op<std::plus<StorageIndex> >(1) 473 //std::bind2nd(std::plus<StorageIndex>(), 1) 474 ); 475 // gsInfo << "ia (outer) :"<< m_matrix.outerSize()<< " / "<< m_matrix.outerSize()+1 <<"\n"; 476 // gsInfo << "ja (inner) :"<< m_matrix.innerSize()<< " / "<< m_matrix.nonZeros() <<"\n"; 477 } 478 }; 479 480 /** \ingroup PardisoSupport_Module 481 * \class PardisoLLT 482 * \brief A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library 483 * 484 * This class allows to solve for A.X = B sparse linear problems via a LL^T Cholesky factorization 485 * using the Intel MKL PARDISO library. The sparse matrix A must be selfajoint and positive definite. 486 * The vectors or matrices X and B can be either dense or sparse. 487 * 488 * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set: 489 * \code solver.pardisoParameterArray()[59] = 1; \endcode 490 * 491 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 492 * \tparam UpLo can be any bitwise combination of Upper, Lower. The default is Upper, meaning only the upper triangular part has to be used. 493 * Upper|Lower can be used to tell both triangular parts can be used as input. 494 * 495 * \implsparsesolverconcept 496 * 497 * \sa \ref TutorialSparseSolverConcept, class SimplicialLLT 498 */ 499 template<typename MatrixType, int _UpLo> 500 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> > 501 { 502 protected: 503 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base; 504 using Base::pardisoInit; 505 using Base::m_matrix; 506 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >; 507 508 public: 509 510 typedef typename Base::Scalar Scalar; 511 typedef typename Base::RealScalar RealScalar; 512 513 typedef typename Base::StorageIndex StorageIndex; 514 enum { UpLo = _UpLo }; 515 using Base::compute; 516 517 PardisoLLT() 518 : Base() 519 { 520 pardisoInit(Base::ScalarIsComplex ? 4 : 2); 521 } 522 523 explicit PardisoLLT(const MatrixType& matrix) 524 : Base() 525 { 526 pardisoInit(Base::ScalarIsComplex ? 4 : 2); 527 compute(matrix); 528 } 529 530 protected: 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 //1-index based conversion 540 StorageIndex * s = m_matrix.innerIndexPtr(); 541 std::transform(s, s+m_matrix.nonZeros(), s, 542 internal::bind2nd_op<std::plus<StorageIndex> >(1) 543 //std::bind2nd(std::plus<StorageIndex>(), 1) 544 ); 545 s = m_matrix.outerIndexPtr(); 546 std::transform(s, s+m_matrix.outerSize()+1, s, 547 internal::bind2nd_op<std::plus<StorageIndex> >(1) 548 //std::bind2nd(std::plus<StorageIndex>(), 1) 549 ); 550 } 551 }; 552 553 /** \ingroup PardisoSupport_Module 554 * \class PardisoLDLT 555 * \brief A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library 556 * 557 * This class allows to solve for A.X = B sparse linear problems via a LDL^T Cholesky factorization 558 * using the Intel MKL PARDISO library. The sparse matrix A is assumed to be selfajoint and positive definite. 559 * For complex matrices, A can also be symmetric only, see the \a Options template parameter. 560 * The vectors or matrices X and B can be either dense or sparse. 561 * 562 * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set: 563 * \code solver.pardisoParameterArray()[59] = 1; \endcode 564 * 565 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 566 * \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. 567 * Symmetric can be used for symmetric, non-selfadjoint complex matrices, the default being to assume a selfadjoint matrix. 568 * Upper|Lower can be used to tell both triangular parts can be used as input. 569 * 570 * \implsparsesolverconcept 571 * 572 * \sa \ref TutorialSparseSolverConcept, class SimplicialLDLT 573 */ 574 template<typename MatrixType, int Options> 575 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > 576 { 577 protected: 578 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base; 579 using Base::pardisoInit; 580 using Base::m_matrix; 581 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >; 582 583 public: 584 585 typedef typename Base::Scalar Scalar; 586 typedef typename Base::RealScalar RealScalar; 587 588 typedef typename Base::StorageIndex StorageIndex; 589 using Base::compute; 590 enum { UpLo = Options&(Upper|Lower) }; 591 592 PardisoLDLT() 593 : Base() 594 { 595 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 596 } 597 598 explicit PardisoLDLT(const MatrixType& matrix) 599 : Base() 600 { 601 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 602 compute(matrix); 603 } 604 605 void getMatrix(const MatrixType& matrix) 606 { 607 // PARDISO supports only upper, row-major matrices 608 //PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null; 609 m_matrix.resize(matrix.rows(), matrix.cols()); 610 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>(); //.twistedBy(p_null); 611 m_matrix.makeCompressed(); 612 //1-index based conversion 613 StorageIndex * s = m_matrix.innerIndexPtr(); 614 std::transform(s, s + m_matrix.nonZeros(), s, 615 internal::bind2nd_op<std::plus<StorageIndex> >(1) 616 //std::bind2nd(std::plus<StorageIndex>(), 1) 617 ); 618 s = m_matrix.outerIndexPtr(); 619 std::transform(s, s + m_matrix.outerSize()+1, s, 620 internal::bind2nd_op<std::plus<StorageIndex> >(1) 621 //std::bind2nd(std::plus<StorageIndex>(), 1) 622 ); 623 } 624 }; 625 626 } // end namespace Eigen 627 628 #endif // EIGEN_PARDISOSUPPORT_H 629