![]() |
Eigen
3.3.3
|
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 IndexType> 00044 struct pardiso_run_selector 00045 { 00046 static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, 00047 IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x) 00048 { 00049 IndexType 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 IndexType; 00058 static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, 00059 IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x) 00060 { 00061 IndexType 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::StorageIndex StorageIndex; 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::StorageIndex StorageIndex; 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::StorageIndex StorageIndex; 00094 }; 00095 00096 } // end namespace internal 00097 00098 template<class Derived> 00099 class PardisoImpl : public SparseSolverBase<Derived> 00100 { 00101 protected: 00102 typedef SparseSolverBase<Derived> Base; 00103 using Base::derived; 00104 using Base::m_isInitialized; 00105 00106 typedef internal::pardiso_traits<Derived> Traits; 00107 public: 00108 using Base::_solve_impl; 00109 00110 typedef typename Traits::MatrixType MatrixType; 00111 typedef typename Traits::Scalar Scalar; 00112 typedef typename Traits::RealScalar RealScalar; 00113 typedef typename Traits::StorageIndex StorageIndex; 00114 typedef SparseMatrix<Scalar,RowMajor,StorageIndex> SparseMatrixType; 00115 typedef Matrix<Scalar,Dynamic,1> VectorType; 00116 typedef Matrix<StorageIndex, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; 00117 typedef Matrix<StorageIndex, MatrixType::RowsAtCompileTime, 1> IntColVectorType; 00118 typedef Array<StorageIndex,64,1,DontAlign> ParameterType; 00119 enum { 00120 ScalarIsComplex = NumTraits<Scalar>::IsComplex, 00121 ColsAtCompileTime = Dynamic, 00122 MaxColsAtCompileTime = Dynamic 00123 }; 00124 00125 PardisoImpl() 00126 { 00127 eigen_assert((sizeof(StorageIndex) >= sizeof(_INTEGER_t) && sizeof(StorageIndex) <= 8) && "Non-supported index type"); 00128 m_iparm.setZero(); 00129 m_msglvl = 0; // No output 00130 m_isInitialized = false; 00131 } 00132 00133 ~PardisoImpl() 00134 { 00135 pardisoRelease(); 00136 } 00137 00138 inline Index cols() const { return m_size; } 00139 inline Index rows() const { return m_size; } 00140 00146 ComputationInfo info() const 00147 { 00148 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00149 return m_info; 00150 } 00151 00155 ParameterType& pardisoParameterArray() 00156 { 00157 return m_iparm; 00158 } 00159 00166 Derived& analyzePattern(const MatrixType& matrix); 00167 00174 Derived& factorize(const MatrixType& matrix); 00175 00176 Derived& compute(const MatrixType& matrix); 00177 00178 template<typename Rhs,typename Dest> 00179 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; 00180 00181 protected: 00182 void pardisoRelease() 00183 { 00184 if(m_isInitialized) // Factorization ran at least once 00185 { 00186 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, 00187 m_iparm.data(), m_msglvl, NULL, NULL); 00188 m_isInitialized = false; 00189 } 00190 } 00191 00192 void pardisoInit(int type) 00193 { 00194 m_type = type; 00195 bool symmetric = std::abs(m_type) < 10; 00196 m_iparm[0] = 1; // No solver default 00197 m_iparm[1] = 2; // use Metis for the ordering 00198 m_iparm[2] = 0; // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??) 00199 m_iparm[3] = 0; // No iterative-direct algorithm 00200 m_iparm[4] = 0; // No user fill-in reducing permutation 00201 m_iparm[5] = 0; // Write solution into x, b is left unchanged 00202 m_iparm[6] = 0; // Not in use 00203 m_iparm[7] = 2; // Max numbers of iterative refinement steps 00204 m_iparm[8] = 0; // Not in use 00205 m_iparm[9] = 13; // Perturb the pivot elements with 1E-13 00206 m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS 00207 m_iparm[11] = 0; // Not in use 00208 m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric). 00209 // Try m_iparm[12] = 1 in case of inappropriate accuracy 00210 m_iparm[13] = 0; // Output: Number of perturbed pivots 00211 m_iparm[14] = 0; // Not in use 00212 m_iparm[15] = 0; // Not in use 00213 m_iparm[16] = 0; // Not in use 00214 m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU 00215 m_iparm[18] = -1; // Output: Mflops for LU factorization 00216 m_iparm[19] = 0; // Output: Numbers of CG Iterations 00217 00218 m_iparm[20] = 0; // 1x1 pivoting 00219 m_iparm[26] = 0; // No matrix checker 00220 m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0; 00221 m_iparm[34] = 1; // C indexing 00222 m_iparm[36] = 0; // CSR 00223 m_iparm[59] = 0; // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core 00224 00225 memset(m_pt, 0, sizeof(m_pt)); 00226 } 00227 00228 protected: 00229 // cached data to reduce reallocation, etc. 00230 00231 void manageErrorCode(Index error) const 00232 { 00233 switch(error) 00234 { 00235 case 0: 00236 m_info = Success; 00237 break; 00238 case -4: 00239 case -7: 00240 m_info = NumericalIssue; 00241 break; 00242 default: 00243 m_info = InvalidInput; 00244 } 00245 } 00246 00247 mutable SparseMatrixType m_matrix; 00248 mutable ComputationInfo m_info; 00249 bool m_analysisIsOk, m_factorizationIsOk; 00250 StorageIndex m_type, m_msglvl; 00251 mutable void *m_pt[64]; 00252 mutable ParameterType m_iparm; 00253 mutable IntColVectorType m_perm; 00254 Index m_size; 00255 00256 }; 00257 00258 template<class Derived> 00259 Derived& PardisoImpl<Derived>::compute(const MatrixType& a) 00260 { 00261 m_size = a.rows(); 00262 eigen_assert(a.rows() == a.cols()); 00263 00264 pardisoRelease(); 00265 m_perm.setZero(m_size); 00266 derived().getMatrix(a); 00267 00268 Index error; 00269 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size), 00270 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 00271 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 00272 manageErrorCode(error); 00273 m_analysisIsOk = true; 00274 m_factorizationIsOk = true; 00275 m_isInitialized = true; 00276 return derived(); 00277 } 00278 00279 template<class Derived> 00280 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a) 00281 { 00282 m_size = a.rows(); 00283 eigen_assert(m_size == a.cols()); 00284 00285 pardisoRelease(); 00286 m_perm.setZero(m_size); 00287 derived().getMatrix(a); 00288 00289 Index error; 00290 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, internal::convert_index<StorageIndex>(m_size), 00291 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 00292 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 00293 00294 manageErrorCode(error); 00295 m_analysisIsOk = true; 00296 m_factorizationIsOk = false; 00297 m_isInitialized = true; 00298 return derived(); 00299 } 00300 00301 template<class Derived> 00302 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a) 00303 { 00304 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 00305 eigen_assert(m_size == a.rows() && m_size == a.cols()); 00306 00307 derived().getMatrix(a); 00308 00309 Index error; 00310 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, internal::convert_index<StorageIndex>(m_size), 00311 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 00312 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL); 00313 00314 manageErrorCode(error); 00315 m_factorizationIsOk = true; 00316 return derived(); 00317 } 00318 00319 template<class Derived> 00320 template<typename BDerived,typename XDerived> 00321 void PardisoImpl<Derived>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const 00322 { 00323 if(m_iparm[0] == 0) // Factorization was not computed 00324 { 00325 m_info = InvalidInput; 00326 return; 00327 } 00328 00329 //Index n = m_matrix.rows(); 00330 Index nrhs = Index(b.cols()); 00331 eigen_assert(m_size==b.rows()); 00332 eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported"); 00333 eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported"); 00334 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows())); 00335 00336 00337 // switch (transposed) { 00338 // case SvNoTrans : m_iparm[11] = 0 ; break; 00339 // case SvTranspose : m_iparm[11] = 2 ; break; 00340 // case SvAdjoint : m_iparm[11] = 1 ; break; 00341 // default: 00342 // //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n"; 00343 // m_iparm[11] = 0; 00344 // } 00345 00346 Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data()); 00347 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp; 00348 00349 // Pardiso cannot solve in-place 00350 if(rhs_ptr == x.derived().data()) 00351 { 00352 tmp = b; 00353 rhs_ptr = tmp.data(); 00354 } 00355 00356 Index error; 00357 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, internal::convert_index<StorageIndex>(m_size), 00358 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(), 00359 m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl, 00360 rhs_ptr, x.derived().data()); 00361 00362 manageErrorCode(error); 00363 } 00364 00365 00383 template<typename MatrixType> 00384 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> > 00385 { 00386 protected: 00387 typedef PardisoImpl<PardisoLU> Base; 00388 typedef typename Base::Scalar Scalar; 00389 typedef typename Base::RealScalar RealScalar; 00390 using Base::pardisoInit; 00391 using Base::m_matrix; 00392 friend class PardisoImpl< PardisoLU<MatrixType> >; 00393 00394 public: 00395 00396 using Base::compute; 00397 using Base::solve; 00398 00399 PardisoLU() 00400 : Base() 00401 { 00402 pardisoInit(Base::ScalarIsComplex ? 13 : 11); 00403 } 00404 00405 explicit PardisoLU(const MatrixType& matrix) 00406 : Base() 00407 { 00408 pardisoInit(Base::ScalarIsComplex ? 13 : 11); 00409 compute(matrix); 00410 } 00411 protected: 00412 void getMatrix(const MatrixType& matrix) 00413 { 00414 m_matrix = matrix; 00415 m_matrix.makeCompressed(); 00416 } 00417 }; 00418 00438 template<typename MatrixType, int _UpLo> 00439 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> > 00440 { 00441 protected: 00442 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base; 00443 typedef typename Base::Scalar Scalar; 00444 typedef typename Base::RealScalar RealScalar; 00445 using Base::pardisoInit; 00446 using Base::m_matrix; 00447 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >; 00448 00449 public: 00450 00451 typedef typename Base::StorageIndex StorageIndex; 00452 enum { UpLo = _UpLo }; 00453 using Base::compute; 00454 00455 PardisoLLT() 00456 : Base() 00457 { 00458 pardisoInit(Base::ScalarIsComplex ? 4 : 2); 00459 } 00460 00461 explicit PardisoLLT(const MatrixType& matrix) 00462 : Base() 00463 { 00464 pardisoInit(Base::ScalarIsComplex ? 4 : 2); 00465 compute(matrix); 00466 } 00467 00468 protected: 00469 00470 void getMatrix(const MatrixType& matrix) 00471 { 00472 // PARDISO supports only upper, row-major matrices 00473 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null; 00474 m_matrix.resize(matrix.rows(), matrix.cols()); 00475 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); 00476 m_matrix.makeCompressed(); 00477 } 00478 }; 00479 00501 template<typename MatrixType, int Options> 00502 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > 00503 { 00504 protected: 00505 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base; 00506 typedef typename Base::Scalar Scalar; 00507 typedef typename Base::RealScalar RealScalar; 00508 using Base::pardisoInit; 00509 using Base::m_matrix; 00510 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >; 00511 00512 public: 00513 00514 typedef typename Base::StorageIndex StorageIndex; 00515 using Base::compute; 00516 enum { UpLo = Options&(Upper|Lower) }; 00517 00518 PardisoLDLT() 00519 : Base() 00520 { 00521 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 00522 } 00523 00524 explicit PardisoLDLT(const MatrixType& matrix) 00525 : Base() 00526 { 00527 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); 00528 compute(matrix); 00529 } 00530 00531 void getMatrix(const MatrixType& matrix) 00532 { 00533 // PARDISO supports only upper, row-major matrices 00534 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null; 00535 m_matrix.resize(matrix.rows(), matrix.cols()); 00536 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null); 00537 m_matrix.makeCompressed(); 00538 } 00539 }; 00540 00541 } // end namespace Eigen 00542 00543 #endif // EIGEN_PARDISOSUPPORT_H