$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 /* 00002 Copyright (c) 2011, Intel Corporation. All rights reserved. 00003 00004 Redistribution and use in source and binary forms, with or without modification, 00005 are permitted provided that the following conditions are met: 00006 00007 * Redistributions of source code must retain the above copyright notice, this 00008 list of conditions and the following disclaimer. 00009 * Redistributions in binary form must reproduce the above copyright notice, 00010 this list of conditions and the following disclaimer in the documentation 00011 and/or other materials provided with the distribution. 00012 * Neither the name of Intel Corporation nor the names of its contributors may 00013 be used to endorse or promote products derived from this software without 00014 specific prior written permission. 00015 00016 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 00017 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 00018 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 00019 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 00020 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 00021 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 00022 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 00023 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00024 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00025 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00026 00027 ******************************************************************************** 00028 * Content : Eigen bindings to Intel(R) MKL PARDISO 00029 ******************************************************************************** 00030 */ 00031 00032 #ifndef EIGEN_PARDISOSUPPORT_H 00033 #define EIGEN_PARDISOSUPPORT_H 00034 00035 namespace Eigen { 00036 00037 template<typename _MatrixType> class PardisoLU; 00038 template<typename _MatrixType, int Options=Upper> class PardisoLLT; 00039 template<typename _MatrixType, int Options=Upper> class PardisoLDLT; 00040 00041 namespace internal 00042 { 00043 template<typename Index> 00044 struct pardiso_run_selector 00045 { 00046 static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a, 00047 Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x) 00048 { 00049 Index error = 0; 00050 ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error); 00051 return error; 00052 } 00053 }; 00054 template<> 00055 struct pardiso_run_selector<long long int> 00056 { 00057 typedef long long int Index; 00058 static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a, 00059 Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x) 00060 { 00061 Index error = 0; 00062 ::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error); 00063 return error; 00064 } 00065 }; 00066 00067 template<class Pardiso> struct pardiso_traits; 00068 00069 template<typename _MatrixType> 00070 struct pardiso_traits< PardisoLU<_MatrixType> > 00071 { 00072 typedef _MatrixType MatrixType; 00073 typedef typename _MatrixType::Scalar Scalar; 00074 typedef typename _MatrixType::RealScalar RealScalar; 00075 typedef typename _MatrixType::Index Index; 00076 }; 00077 00078 template<typename _MatrixType, int Options> 00079 struct pardiso_traits< PardisoLLT<_MatrixType, Options> > 00080 { 00081 typedef _MatrixType MatrixType; 00082 typedef typename _MatrixType::Scalar Scalar; 00083 typedef typename _MatrixType::RealScalar RealScalar; 00084 typedef typename _MatrixType::Index Index; 00085 }; 00086 00087 template<typename _MatrixType, int Options> 00088 struct pardiso_traits< PardisoLDLT<_MatrixType, Options> > 00089 { 00090 typedef _MatrixType MatrixType; 00091 typedef typename _MatrixType::Scalar Scalar; 00092 typedef typename _MatrixType::RealScalar RealScalar; 00093 typedef typename _MatrixType::Index Index; 00094 }; 00095 00096 } 00097 00098 template<class Derived> 00099 class PardisoImpl 00100 { 00101 typedef internal::pardiso_traits<Derived> Traits; 00102 public: 00103 typedef typename Traits::MatrixType MatrixType; 00104 typedef typename Traits::Scalar Scalar; 00105 typedef typename Traits::RealScalar RealScalar; 00106 typedef typename Traits::Index Index; 00107 typedef SparseMatrix<Scalar,RowMajor,Index> SparseMatrixType; 00108 typedef Matrix<Scalar,Dynamic,1> VectorType; 00109 typedef Matrix<Index, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 00110 typedef Matrix<Index, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 00111 typedef Array<Index,64,1,DontAlign> ParameterType; 00112 enum { 00113 ScalarIsComplex = NumTraits<Scalar>::IsComplex 00114 }; 00115 00116 PardisoImpl() 00117 { 00118 eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type"); 00119 m_iparm.setZero(); 00120 m_msglvl = 0; // No output 00121 m_initialized = false; 00122 } 00123 00124 ~PardisoImpl() 00125 { 00126 pardisoRelease(); 00127 } 00128 00129 inline Index cols() const { return m_size; } 00130 inline Index rows() const { return m_size; } 00131 00137 ComputationInfo info() const 00138 { 00139 eigen_assert(m_initialized && "Decomposition is not initialized."); 00140 return m_info; 00141 } 00142 00146 ParameterType& pardisoParameterArray() 00147 { 00148 return m_iparm; 00149 } 00150 00157 Derived& analyzePattern(const MatrixType& matrix); 00158 00165 Derived& factorize(const MatrixType& matrix); 00166 00167 Derived& compute(const MatrixType& matrix); 00168 00173 template<typename Rhs> 00174 inline const internal::solve_retval<PardisoImpl, Rhs> 00175 solve(const MatrixBase<Rhs>& b) const 00176 { 00177 eigen_assert(m_initialized && "Pardiso solver is not initialized."); 00178 eigen_assert(rows()==b.rows() 00179 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b"); 00180 return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived()); 00181 } 00182 00187 template<typename Rhs> 00188 inline const internal::sparse_solve_retval<PardisoImpl, Rhs> 00189 solve(const SparseMatrixBase<Rhs>& b) const 00190 { 00191 eigen_assert(m_initialized && "Pardiso solver is not initialized."); 00192 eigen_assert(rows()==b.rows() 00193 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b"); 00194 return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived()); 00195 } 00196 00197 Derived& derived() 00198 { 00199 return *static_cast<Derived*>(this); 00200 } 00201 const Derived& derived() const 00202 { 00203 return *static_cast<const Derived*>(this); 00204 } 00205 00206 template<typename BDerived, typename XDerived> 00207 bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const; 00208 00209 protected: 00210 void pardisoRelease() 00211 { 00212 if(m_initialized) // Factorization ran at least once 00213 { 00214 internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0, 00215 m_iparm.data(), m_msglvl, 0, 0); 00216 } 00217 } 00218 00219 void pardisoInit(int type) 00220 { 00221 m_type = type; 00222 bool symmetric = std::abs(m_type) < 10; 00223 m_iparm[0] = 1; // No solver default 00224 m_iparm[1] = 3; // use Metis for the ordering 00225 m_iparm[2] = 1; // Numbers of processors, value of OMP_NUM_THREADS 00226 m_iparm[3] = 0; // No iterative-direct algorithm 00227 m_iparm[4] = 0; // No user fill-in reducing permutation 00228 m_iparm[5] = 0; // Write solution into x 00229 m_iparm[6] = 0; // Not in use 00230 m_iparm[7] = 2; // Max numbers of iterative refinement steps 00231 m_iparm[8] = 0; // Not in use 00232 m_iparm[9] = 13; // Perturb the pivot elements with 1E-13 00233 m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS 00234 m_iparm[11] = 0; // Not in use 00235 m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric). 00236 // Try m_iparm[12] = 1 in case of inappropriate accuracy 00237 m_iparm[13] = 0; // Output: Number of perturbed pivots 00238 m_iparm[14] = 0; // Not in use 00239 m_iparm[15] = 0; // Not in use 00240 m_iparm[16] = 0; // Not in use 00241 m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU 00242 m_iparm[18] = -1; // Output: Mflops for LU factorization 00243 m_iparm[19] = 0; // Output: Numbers of CG Iterations 00244 00245 m_iparm[20] = 0; // 1x1 pivoting 00246 m_iparm[26] = 0; // No matrix checker 00247 m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0; 00248 m_iparm[34] = 1; // C indexing 00249 m_iparm[59] = 1; // Automatic switch between In-Core and Out-of-Core modes 00250 } 00251 00252 protected: 00253 // cached data to reduce reallocation, etc. 00254 00255 void manageErrorCode(Index error) 00256 { 00257 switch(error) 00258 { 00259 case 0: 00260 m_info = Success; 00261 break; 00262 case -4: 00263 case -7: 00264 m_info = NumericalIssue; 00265 break; 00266 default: 00267 m_info = InvalidInput; 00268 } 00269 } 00270 00271 mutable SparseMatrixType m_matrix; 00272 ComputationInfo m_info; 00273 bool m_initialized, m_analysisIsOk, m_factorizationIsOk; 00274 Index m_type, m_msglvl; 00275 mutable void *m_pt[64]; 00276 mutable ParameterType m_iparm; 00277 mutable IntColVectorType m_perm; 00278 Index m_size; 00279 00280 private: 00281 PardisoImpl(PardisoImpl &) {} 00282 }; 00283 00284 template<class Derived> 00285 Derived& PardisoImpl<Derived>::compute(const MatrixType& a) 00286 { 00287 m_size = a.rows(); 00288 eigen_assert(a.rows() == a.cols()); 00289 00290 pardisoRelease(); 00291 memset(m_pt, 0, sizeof(m_pt)); 00292 m_perm.setZero(m_size); 00293 derived().getMatrix(a); 00294 00295 Index error; 00296 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size, 00297 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 00298 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 00299 00300 manageErrorCode(error); 00301 m_analysisIsOk = true; 00302 m_factorizationIsOk = true; 00303 m_initialized = true; 00304 return derived(); 00305 } 00306 00307 template<class Derived> 00308 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a) 00309 { 00310 m_size = a.rows(); 00311 eigen_assert(m_size == a.cols()); 00312 00313 pardisoRelease(); 00314 memset(m_pt, 0, sizeof(m_pt)); 00315 m_perm.setZero(m_size); 00316 derived().getMatrix(a); 00317 00318 Index error; 00319 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size, 00320 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 00321 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 00322 00323 manageErrorCode(error); 00324 m_analysisIsOk = true; 00325 m_factorizationIsOk = false; 00326 m_initialized = true; 00327 return derived(); 00328 } 00329 00330 template<class Derived> 00331 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a) 00332 { 00333 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 00334 eigen_assert(m_size == a.rows() && m_size == a.cols()); 00335 00336 derived().getMatrix(a); 00337 00338 Index error; 00339 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size, 00340 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 00341 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 00342 00343 manageErrorCode(error); 00344 m_factorizationIsOk = true; 00345 return derived(); 00346 } 00347 00348 template<class Base> 00349 template<typename BDerived,typename XDerived> 00350 bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const 00351 { 00352 if(m_iparm[0] == 0) // Factorization was not computed 00353 return false; 00354 00355 //Index n = m_matrix.rows(); 00356 Index nrhs = Index(b.cols()); 00357 eigen_assert(m_size==b.rows()); 00358 eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported"); 00359 eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported"); 00360 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows())); 00361 00362 00363 // switch (transposed) { 00364 // case SvNoTrans : m_iparm[11] = 0 ; break; 00365 // case SvTranspose : m_iparm[11] = 2 ; break; 00366 // case SvAdjoint : m_iparm[11] = 1 ; break; 00367 // default: 00368 // //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n"; 00369 // m_iparm[11] = 0; 00370 // } 00371 00372 Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data()); 00373 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp; 00374 00375 // Pardiso cannot solve in-place 00376 if(rhs_ptr == x.derived().data()) 00377 { 00378 tmp = b; 00379 rhs_ptr = tmp.data(); 00380 } 00381 00382 Index error; 00383 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size, 00384 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 00385 m_perm.data(), nrhs, m_iparm.data(), m_msglvl, 00386 rhs_ptr, x.derived().data()); 00387 00388 return error==0; 00389 } 00390 00391 00404 template<typename MatrixType> 00405 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> > 00406 { 00407 protected: 00408 typedef PardisoImpl< PardisoLU<MatrixType> > Base; 00409 typedef typename Base::Scalar Scalar; 00410 typedef typename Base::RealScalar RealScalar; 00411 using Base::pardisoInit; 00412 using Base::m_matrix; 00413 friend class PardisoImpl< PardisoLU<MatrixType> >; 00414 00415 public: 00416 00417 using Base::compute; 00418 using Base::solve; 00419 00420 PardisoLU() 00421 : Base() 00422 { 00423 pardisoInit(Base::ScalarIsComplex ? 13 : 11); 00424 } 00425 00426 PardisoLU(const MatrixType& matrix) 00427 : Base() 00428 { 00429 pardisoInit(Base::ScalarIsComplex ? 13 : 11); 00430 compute(matrix); 00431 } 00432 protected: 00433 void getMatrix(const MatrixType& matrix) 00434 { 00435 m_matrix = matrix; 00436 } 00437 00438 private: 00439 PardisoLU(PardisoLU& ) {} 00440 }; 00441 00456 template<typename MatrixType, int _UpLo> 00457 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> > 00458 { 00459 protected: 00460 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base; 00461 typedef typename Base::Scalar Scalar; 00462 typedef typename Base::Index Index; 00463 typedef typename Base::RealScalar RealScalar; 00464 using Base::pardisoInit; 00465 using Base::m_matrix; 00466 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >; 00467 00468 public: 00469 00470 enum { UpLo = _UpLo }; 00471 using Base::compute; 00472 using Base::solve; 00473 00474 PardisoLLT() 00475 : Base() 00476 { 00477 pardisoInit(Base::ScalarIsComplex ? 4 : 2); 00478 } 00479 00480 PardisoLLT(const MatrixType& matrix) 00481 : Base() 00482 { 00483 pardisoInit(Base::ScalarIsComplex ? 4 : 2); 00484 compute(matrix); 00485 } 00486 00487 protected: 00488 00489 void getMatrix(const MatrixType& matrix) 00490 { 00491 // PARDISO supports only upper, row-major matrices 00492 PermutationMatrix<Dynamic,Dynamic,Index> p_null; 00493 m_matrix.resize(matrix.rows(), matrix.cols()); 00494 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); 00495 } 00496 00497 private: 00498 PardisoLLT(PardisoLLT& ) {} 00499 }; 00500 00517 template<typename MatrixType, int Options> 00518 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > 00519 { 00520 protected: 00521 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base; 00522 typedef typename Base::Scalar Scalar; 00523 typedef typename Base::Index Index; 00524 typedef typename Base::RealScalar RealScalar; 00525 using Base::pardisoInit; 00526 using Base::m_matrix; 00527 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >; 00528 00529 public: 00530 00531 using Base::compute; 00532 using Base::solve; 00533 enum { UpLo = Options&(Upper|Lower) }; 00534 00535 PardisoLDLT() 00536 : Base() 00537 { 00538 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 00539 } 00540 00541 PardisoLDLT(const MatrixType& matrix) 00542 : Base() 00543 { 00544 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 00545 compute(matrix); 00546 } 00547 00548 void getMatrix(const MatrixType& matrix) 00549 { 00550 // PARDISO supports only upper, row-major matrices 00551 PermutationMatrix<Dynamic,Dynamic,Index> p_null; 00552 m_matrix.resize(matrix.rows(), matrix.cols()); 00553 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); 00554 } 00555 00556 private: 00557 PardisoLDLT(PardisoLDLT& ) {} 00558 }; 00559 00560 namespace internal { 00561 00562 template<typename _Derived, typename Rhs> 00563 struct solve_retval<PardisoImpl<_Derived>, Rhs> 00564 : solve_retval_base<PardisoImpl<_Derived>, Rhs> 00565 { 00566 typedef PardisoImpl<_Derived> Dec; 00567 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00568 00569 template<typename Dest> void evalTo(Dest& dst) const 00570 { 00571 dec()._solve(rhs(),dst); 00572 } 00573 }; 00574 00575 template<typename Derived, typename Rhs> 00576 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs> 00577 : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs> 00578 { 00579 typedef PardisoImpl<Derived> Dec; 00580 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 00581 00582 template<typename Dest> void evalTo(Dest& dst) const 00583 { 00584 this->defaultEvalTo(dst); 00585 } 00586 }; 00587 00588 } // end namespace internal 00589 00590 } // end namespace Eigen 00591 00592 #endif // EIGEN_PARDISOSUPPORT_H