32 #ifndef EIGEN_PARDISOSUPPORT_H
33 #define EIGEN_PARDISOSUPPORT_H
38 template<
typename _MatrixType,
int Options=Upper>
class PardisoLLT;
39 template<
typename _MatrixType,
int Options=Upper>
class PardisoLDLT;
43 template<
typename IndexType>
44 struct pardiso_run_selector
46 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)
50 ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
55 struct pardiso_run_selector<long long int>
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)
62 ::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
67 template<
class Pardiso>
struct pardiso_traits;
69 template<
typename _MatrixType>
70 struct pardiso_traits< PardisoLU<_MatrixType> >
72 typedef _MatrixType MatrixType;
73 typedef typename _MatrixType::Scalar Scalar;
74 typedef typename _MatrixType::RealScalar RealScalar;
75 typedef typename _MatrixType::StorageIndex StorageIndex;
78 template<
typename _MatrixType,
int Options>
79 struct pardiso_traits< PardisoLLT<_MatrixType, Options> >
81 typedef _MatrixType MatrixType;
82 typedef typename _MatrixType::Scalar Scalar;
83 typedef typename _MatrixType::RealScalar RealScalar;
84 typedef typename _MatrixType::StorageIndex StorageIndex;
87 template<
typename _MatrixType,
int Options>
88 struct pardiso_traits< PardisoLDLT<_MatrixType, Options> >
90 typedef _MatrixType MatrixType;
91 typedef typename _MatrixType::Scalar Scalar;
92 typedef typename _MatrixType::RealScalar RealScalar;
93 typedef typename _MatrixType::StorageIndex StorageIndex;
98 template<
class Derived>
99 class PardisoImpl :
public SparseSolverBase<Derived>
102 typedef SparseSolverBase<Derived> Base;
104 using Base::m_isInitialized;
106 typedef internal::pardiso_traits<Derived> Traits;
108 using Base::_solve_impl;
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;
120 ScalarIsComplex = NumTraits<Scalar>::IsComplex
125 eigen_assert((
sizeof(StorageIndex) >=
sizeof(_INTEGER_t) &&
sizeof(StorageIndex) <= 8) &&
"Non-supported index type");
128 m_isInitialized =
false;
136 inline Index cols()
const {
return m_size; }
137 inline Index rows()
const {
return m_size; }
146 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
153 ParameterType& pardisoParameterArray()
164 Derived& analyzePattern(
const MatrixType& matrix);
172 Derived& factorize(
const MatrixType& matrix);
174 Derived& compute(
const MatrixType& matrix);
176 template<
typename Rhs,
typename Dest>
177 void _solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const;
180 void pardisoRelease()
184 internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, -1, m_size,0, 0, 0, m_perm.data(), 0,
185 m_iparm.data(), m_msglvl, NULL, NULL);
186 m_isInitialized =
false;
190 void pardisoInit(
int type)
193 bool symmetric = std::abs(m_type) < 10;
204 m_iparm[10] = symmetric ? 0 : 1;
206 m_iparm[12] = symmetric ? 0 : 1;
218 m_iparm[27] = (
sizeof(RealScalar) == 4) ? 1 : 0;
222 memset(m_pt, 0,
sizeof(m_pt));
228 void manageErrorCode(Index error)
const
244 mutable SparseMatrixType m_matrix;
246 bool m_analysisIsOk, m_factorizationIsOk;
247 Index m_type, m_msglvl;
248 mutable void *m_pt[64];
249 mutable ParameterType m_iparm;
250 mutable IntColVectorType m_perm;
255 template<
class Derived>
256 Derived& PardisoImpl<Derived>::compute(
const MatrixType& a)
259 eigen_assert(a.rows() == a.cols());
262 m_perm.setZero(m_size);
263 derived().getMatrix(a);
266 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, m_size,
267 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
268 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
270 manageErrorCode(error);
271 m_analysisIsOk =
true;
272 m_factorizationIsOk =
true;
273 m_isInitialized =
true;
277 template<
class Derived>
278 Derived& PardisoImpl<Derived>::analyzePattern(
const MatrixType& a)
281 eigen_assert(m_size == a.cols());
284 m_perm.setZero(m_size);
285 derived().getMatrix(a);
288 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, m_size,
289 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
290 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
292 manageErrorCode(error);
293 m_analysisIsOk =
true;
294 m_factorizationIsOk =
false;
295 m_isInitialized =
true;
299 template<
class Derived>
300 Derived& PardisoImpl<Derived>::factorize(
const MatrixType& a)
302 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
303 eigen_assert(m_size == a.rows() && m_size == a.cols());
305 derived().getMatrix(a);
308 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, m_size,
309 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
310 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
312 manageErrorCode(error);
313 m_factorizationIsOk =
true;
317 template<
class Derived>
318 template<
typename BDerived,
typename XDerived>
319 void PardisoImpl<Derived>::_solve_impl(
const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x)
const
328 Index nrhs = Index(b.cols());
329 eigen_assert(m_size==b.rows());
330 eigen_assert(((MatrixBase<BDerived>::Flags &
RowMajorBit) == 0 || nrhs == 1) &&
"Row-major right hand sides are not supported");
331 eigen_assert(((MatrixBase<XDerived>::Flags &
RowMajorBit) == 0 || nrhs == 1) &&
"Row-major matrices of unknowns are not supported");
332 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
344 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.derived().data());
345 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
348 if(rhs_ptr == x.derived().data())
351 rhs_ptr = tmp.data();
355 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, m_size,
356 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
357 m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
358 rhs_ptr, x.derived().data());
360 manageErrorCode(error);
376 template<
typename MatrixType>
377 class PardisoLU :
public PardisoImpl< PardisoLU<MatrixType> >
380 typedef PardisoImpl<PardisoLU> Base;
381 typedef typename Base::Scalar Scalar;
382 typedef typename Base::RealScalar RealScalar;
383 using Base::pardisoInit;
384 using Base::m_matrix;
385 friend class PardisoImpl< PardisoLU<MatrixType> >;
395 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
398 explicit PardisoLU(
const MatrixType& matrix)
401 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
405 void getMatrix(
const MatrixType& matrix)
408 m_matrix.makeCompressed();
426 template<
typename MatrixType,
int _UpLo>
427 class PardisoLLT :
public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
430 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
431 typedef typename Base::Scalar Scalar;
432 typedef typename Base::RealScalar RealScalar;
433 using Base::pardisoInit;
434 using Base::m_matrix;
435 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
439 typedef typename Base::StorageIndex StorageIndex;
440 enum { UpLo = _UpLo };
446 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
449 explicit PardisoLLT(
const MatrixType& matrix)
452 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
458 void getMatrix(
const MatrixType& matrix)
461 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
462 m_matrix.resize(matrix.rows(), matrix.cols());
463 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
464 m_matrix.makeCompressed();
484 template<
typename MatrixType,
int Options>
485 class PardisoLDLT :
public PardisoImpl< PardisoLDLT<MatrixType,Options> >
488 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
489 typedef typename Base::Scalar Scalar;
490 typedef typename Base::RealScalar RealScalar;
491 using Base::pardisoInit;
492 using Base::m_matrix;
493 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
497 typedef typename Base::StorageIndex StorageIndex;
504 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
507 explicit PardisoLDLT(
const MatrixType& matrix)
510 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
514 void getMatrix(
const MatrixType& matrix)
517 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
518 m_matrix.resize(matrix.rows(), matrix.cols());
519 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
520 m_matrix.makeCompressed();
526 #endif // EIGEN_PARDISOSUPPORT_H
Definition: Constants.h:196
Definition: Constants.h:214
const unsigned int RowMajorBit
Definition: Constants.h:53
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:38
Definition: Constants.h:198
Definition: Constants.h:426
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
Definition: Constants.h:431
Definition: Constants.h:424
Definition: Eigen_Colamd.h:54
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:39
ComputationInfo
Definition: Constants.h:422