00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
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 enum {
00112 ScalarIsComplex = NumTraits<Scalar>::IsComplex
00113 };
00114
00115 PardisoImpl()
00116 {
00117 eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type");
00118 m_iparm.setZero();
00119 m_msglvl = 0;
00120 m_initialized = false;
00121 }
00122
00123 ~PardisoImpl()
00124 {
00125 pardisoRelease();
00126 }
00127
00128 inline Index cols() const { return m_size; }
00129 inline Index rows() const { return m_size; }
00130
00136 ComputationInfo info() const
00137 {
00138 eigen_assert(m_initialized && "Decomposition is not initialized.");
00139 return m_info;
00140 }
00141
00145 Array<Index,64,1>& pardisoParameterArray()
00146 {
00147 return m_iparm;
00148 }
00149
00156 Derived& analyzePattern(const MatrixType& matrix);
00157
00164 Derived& factorize(const MatrixType& matrix);
00165
00166 Derived& compute(const MatrixType& matrix);
00167
00172 template<typename Rhs>
00173 inline const internal::solve_retval<PardisoImpl, Rhs>
00174 solve(const MatrixBase<Rhs>& b) const
00175 {
00176 eigen_assert(m_initialized && "Pardiso solver is not initialized.");
00177 eigen_assert(rows()==b.rows()
00178 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
00179 return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived());
00180 }
00181
00186 template<typename Rhs>
00187 inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
00188 solve(const SparseMatrixBase<Rhs>& b) const
00189 {
00190 eigen_assert(m_initialized && "Pardiso solver is not initialized.");
00191 eigen_assert(rows()==b.rows()
00192 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
00193 return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived());
00194 }
00195
00196 Derived& derived()
00197 {
00198 return *static_cast<Derived*>(this);
00199 }
00200 const Derived& derived() const
00201 {
00202 return *static_cast<const Derived*>(this);
00203 }
00204
00205 template<typename BDerived, typename XDerived>
00206 bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const;
00207
00209 template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00210 void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00211 {
00212 eigen_assert(m_size==b.rows());
00213
00214
00215 static const int NbColsAtOnce = 4;
00216 int rhsCols = b.cols();
00217 int size = b.rows();
00218
00219
00220 Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_rhs(size,rhsCols);
00221 Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_res(size,rhsCols);
00222 for(int k=0; k<rhsCols; k+=NbColsAtOnce)
00223 {
00224 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
00225 tmp_rhs.leftCols(actualCols) = b.middleCols(k,actualCols);
00226 tmp_res.leftCols(actualCols) = derived().solve(tmp_rhs.leftCols(actualCols));
00227 dest.middleCols(k,actualCols) = tmp_res.leftCols(actualCols).sparseView();
00228 }
00229 }
00230
00231 protected:
00232 void pardisoRelease()
00233 {
00234 if(m_initialized)
00235 {
00236 internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
00237 m_iparm.data(), m_msglvl, 0, 0);
00238 }
00239 }
00240
00241 void pardisoInit(int type)
00242 {
00243 m_type = type;
00244 bool symmetric = abs(m_type) < 10;
00245 m_iparm[0] = 1;
00246 m_iparm[1] = 3;
00247 m_iparm[2] = 1;
00248 m_iparm[3] = 0;
00249 m_iparm[4] = 0;
00250 m_iparm[5] = 0;
00251 m_iparm[6] = 0;
00252 m_iparm[7] = 2;
00253 m_iparm[8] = 0;
00254 m_iparm[9] = 13;
00255 m_iparm[10] = symmetric ? 0 : 1;
00256 m_iparm[11] = 0;
00257 m_iparm[12] = symmetric ? 0 : 1;
00258
00259 m_iparm[13] = 0;
00260 m_iparm[14] = 0;
00261 m_iparm[15] = 0;
00262 m_iparm[16] = 0;
00263 m_iparm[17] = -1;
00264 m_iparm[18] = -1;
00265 m_iparm[19] = 0;
00266
00267 m_iparm[20] = 0;
00268 m_iparm[26] = 0;
00269 m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
00270 m_iparm[34] = 1;
00271 m_iparm[59] = 1;
00272 }
00273
00274 protected:
00275
00276
00277 void manageErrorCode(Index error)
00278 {
00279 switch(error)
00280 {
00281 case 0:
00282 m_info = Success;
00283 break;
00284 case -4:
00285 case -7:
00286 m_info = NumericalIssue;
00287 break;
00288 default:
00289 m_info = InvalidInput;
00290 }
00291 }
00292
00293 mutable SparseMatrixType m_matrix;
00294 ComputationInfo m_info;
00295 bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
00296 Index m_type, m_msglvl;
00297 mutable void *m_pt[64];
00298 mutable Array<Index,64,1> m_iparm;
00299 mutable IntColVectorType m_perm;
00300 Index m_size;
00301
00302 private:
00303 PardisoImpl(PardisoImpl &) {}
00304 };
00305
00306 template<class Derived>
00307 Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
00308 {
00309 m_size = a.rows();
00310 eigen_assert(a.rows() == a.cols());
00311
00312 pardisoRelease();
00313 memset(m_pt, 0, sizeof(m_pt));
00314 m_perm.setZero(m_size);
00315 derived().getMatrix(a);
00316
00317 Index error;
00318 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
00319 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00320 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00321
00322 manageErrorCode(error);
00323 m_analysisIsOk = true;
00324 m_factorizationIsOk = true;
00325 m_initialized = true;
00326 return derived();
00327 }
00328
00329 template<class Derived>
00330 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
00331 {
00332 m_size = a.rows();
00333 eigen_assert(m_size == a.cols());
00334
00335 pardisoRelease();
00336 memset(m_pt, 0, sizeof(m_pt));
00337 m_perm.setZero(m_size);
00338 derived().getMatrix(a);
00339
00340 Index error;
00341 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
00342 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00343 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00344
00345 manageErrorCode(error);
00346 m_analysisIsOk = true;
00347 m_factorizationIsOk = false;
00348 m_initialized = true;
00349 return derived();
00350 }
00351
00352 template<class Derived>
00353 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
00354 {
00355 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00356 eigen_assert(m_size == a.rows() && m_size == a.cols());
00357
00358 derived().getMatrix(a);
00359
00360 Index error;
00361 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
00362 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00363 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00364
00365 manageErrorCode(error);
00366 m_factorizationIsOk = true;
00367 return derived();
00368 }
00369
00370 template<class Base>
00371 template<typename BDerived,typename XDerived>
00372 bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
00373 {
00374 if(m_iparm[0] == 0)
00375 return false;
00376
00377
00378 Index nrhs = Index(b.cols());
00379 eigen_assert(m_size==b.rows());
00380 eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
00381 eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
00382 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394 Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
00395 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
00396
00397
00398 if(rhs_ptr == x.derived().data())
00399 {
00400 tmp = b;
00401 rhs_ptr = tmp.data();
00402 }
00403
00404 Index error;
00405 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
00406 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00407 m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
00408 rhs_ptr, x.derived().data());
00409
00410 return error==0;
00411 }
00412
00413
00426 template<typename MatrixType>
00427 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
00428 {
00429 protected:
00430 typedef PardisoImpl< PardisoLU<MatrixType> > Base;
00431 typedef typename Base::Scalar Scalar;
00432 typedef typename Base::RealScalar RealScalar;
00433 using Base::pardisoInit;
00434 using Base::m_matrix;
00435 friend class PardisoImpl< PardisoLU<MatrixType> >;
00436
00437 public:
00438
00439 using Base::compute;
00440 using Base::solve;
00441
00442 PardisoLU()
00443 : Base()
00444 {
00445 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00446 }
00447
00448 PardisoLU(const MatrixType& matrix)
00449 : Base()
00450 {
00451 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00452 compute(matrix);
00453 }
00454 protected:
00455 void getMatrix(const MatrixType& matrix)
00456 {
00457 m_matrix = matrix;
00458 }
00459
00460 private:
00461 PardisoLU(PardisoLU& ) {}
00462 };
00463
00478 template<typename MatrixType, int _UpLo>
00479 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
00480 {
00481 protected:
00482 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
00483 typedef typename Base::Scalar Scalar;
00484 typedef typename Base::Index Index;
00485 typedef typename Base::RealScalar RealScalar;
00486 using Base::pardisoInit;
00487 using Base::m_matrix;
00488 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
00489
00490 public:
00491
00492 enum { UpLo = _UpLo };
00493 using Base::compute;
00494 using Base::solve;
00495
00496 PardisoLLT()
00497 : Base()
00498 {
00499 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00500 }
00501
00502 PardisoLLT(const MatrixType& matrix)
00503 : Base()
00504 {
00505 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00506 compute(matrix);
00507 }
00508
00509 protected:
00510
00511 void getMatrix(const MatrixType& matrix)
00512 {
00513
00514 PermutationMatrix<Dynamic,Dynamic,Index> p_null;
00515 m_matrix.resize(matrix.rows(), matrix.cols());
00516 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00517 }
00518
00519 private:
00520 PardisoLLT(PardisoLLT& ) {}
00521 };
00522
00539 template<typename MatrixType, int Options>
00540 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
00541 {
00542 protected:
00543 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
00544 typedef typename Base::Scalar Scalar;
00545 typedef typename Base::Index Index;
00546 typedef typename Base::RealScalar RealScalar;
00547 using Base::pardisoInit;
00548 using Base::m_matrix;
00549 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
00550
00551 public:
00552
00553 using Base::compute;
00554 using Base::solve;
00555 enum { UpLo = Options&(Upper|Lower) };
00556
00557 PardisoLDLT()
00558 : Base()
00559 {
00560 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00561 }
00562
00563 PardisoLDLT(const MatrixType& matrix)
00564 : Base()
00565 {
00566 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00567 compute(matrix);
00568 }
00569
00570 void getMatrix(const MatrixType& matrix)
00571 {
00572
00573 PermutationMatrix<Dynamic,Dynamic,Index> p_null;
00574 m_matrix.resize(matrix.rows(), matrix.cols());
00575 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00576 }
00577
00578 private:
00579 PardisoLDLT(PardisoLDLT& ) {}
00580 };
00581
00582 namespace internal {
00583
00584 template<typename _Derived, typename Rhs>
00585 struct solve_retval<PardisoImpl<_Derived>, Rhs>
00586 : solve_retval_base<PardisoImpl<_Derived>, Rhs>
00587 {
00588 typedef PardisoImpl<_Derived> Dec;
00589 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00590
00591 template<typename Dest> void evalTo(Dest& dst) const
00592 {
00593 dec()._solve(rhs(),dst);
00594 }
00595 };
00596
00597 template<typename Derived, typename Rhs>
00598 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
00599 : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
00600 {
00601 typedef PardisoImpl<Derived> Dec;
00602 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00603
00604 template<typename Dest> void evalTo(Dest& dst) const
00605 {
00606 dec().derived()._solve_sparse(rhs(),dst);
00607 }
00608 };
00609
00610 }
00611
00612 }
00613
00614 #endif // EIGEN_PARDISOSUPPORT_H