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,
121 ColsAtCompileTime = Dynamic,
122 MaxColsAtCompileTime = Dynamic
127 eigen_assert((
sizeof(StorageIndex) >=
sizeof(_INTEGER_t) &&
sizeof(StorageIndex) <= 8) &&
"Non-supported index type");
130 m_isInitialized =
false;
138 inline Index cols()
const {
return m_size; }
139 inline Index rows()
const {
return m_size; }
148 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
155 ParameterType& pardisoParameterArray()
166 Derived& analyzePattern(
const MatrixType& matrix);
174 Derived& factorize(
const MatrixType& matrix);
176 Derived& compute(
const MatrixType& matrix);
178 template<
typename Rhs,
typename Dest>
179 void _solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const;
182 void pardisoRelease()
186 internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, -1, m_size,0, 0, 0, m_perm.data(), 0,
187 m_iparm.data(), m_msglvl, NULL, NULL);
188 m_isInitialized =
false;
192 void pardisoInit(
int type)
195 bool symmetric = std::abs(m_type) < 10;
206 m_iparm[10] = symmetric ? 0 : 1;
208 m_iparm[12] = symmetric ? 0 : 1;
220 m_iparm[27] = (
sizeof(RealScalar) == 4) ? 1 : 0;
224 memset(m_pt, 0,
sizeof(m_pt));
230 void manageErrorCode(Index error)
const
246 mutable SparseMatrixType m_matrix;
248 bool m_analysisIsOk, m_factorizationIsOk;
249 Index m_type, m_msglvl;
250 mutable void *m_pt[64];
251 mutable ParameterType m_iparm;
252 mutable IntColVectorType m_perm;
257 template<
class Derived>
258 Derived& PardisoImpl<Derived>::compute(
const MatrixType& a)
261 eigen_assert(a.rows() == a.cols());
264 m_perm.setZero(m_size);
265 derived().getMatrix(a);
268 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, m_size,
269 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
270 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
272 manageErrorCode(error);
273 m_analysisIsOk =
true;
274 m_factorizationIsOk =
true;
275 m_isInitialized =
true;
279 template<
class Derived>
280 Derived& PardisoImpl<Derived>::analyzePattern(
const MatrixType& a)
283 eigen_assert(m_size == a.cols());
286 m_perm.setZero(m_size);
287 derived().getMatrix(a);
290 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, m_size,
291 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
292 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
294 manageErrorCode(error);
295 m_analysisIsOk =
true;
296 m_factorizationIsOk =
false;
297 m_isInitialized =
true;
301 template<
class Derived>
302 Derived& PardisoImpl<Derived>::factorize(
const MatrixType& a)
304 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
305 eigen_assert(m_size == a.rows() && m_size == a.cols());
307 derived().getMatrix(a);
310 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, m_size,
311 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
312 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
314 manageErrorCode(error);
315 m_factorizationIsOk =
true;
319 template<
class Derived>
320 template<
typename BDerived,
typename XDerived>
321 void PardisoImpl<Derived>::_solve_impl(
const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x)
const
330 Index nrhs = Index(b.cols());
331 eigen_assert(m_size==b.rows());
332 eigen_assert(((MatrixBase<BDerived>::Flags &
RowMajorBit) == 0 || nrhs == 1) &&
"Row-major right hand sides are not supported");
333 eigen_assert(((MatrixBase<XDerived>::Flags &
RowMajorBit) == 0 || nrhs == 1) &&
"Row-major matrices of unknowns are not supported");
334 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
346 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.derived().data());
347 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
350 if(rhs_ptr == x.derived().data())
353 rhs_ptr = tmp.data();
357 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, m_size,
358 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
359 m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
360 rhs_ptr, x.derived().data());
362 manageErrorCode(error);
380 template<
typename MatrixType>
381 class PardisoLU :
public PardisoImpl< PardisoLU<MatrixType> >
384 typedef PardisoImpl<PardisoLU> Base;
385 typedef typename Base::Scalar Scalar;
386 typedef typename Base::RealScalar RealScalar;
387 using Base::pardisoInit;
388 using Base::m_matrix;
389 friend class PardisoImpl< PardisoLU<MatrixType> >;
399 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
402 explicit PardisoLU(
const MatrixType& matrix)
405 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
409 void getMatrix(
const MatrixType& matrix)
412 m_matrix.makeCompressed();
432 template<
typename MatrixType,
int _UpLo>
433 class PardisoLLT :
public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
436 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
437 typedef typename Base::Scalar Scalar;
438 typedef typename Base::RealScalar RealScalar;
439 using Base::pardisoInit;
440 using Base::m_matrix;
441 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
445 typedef typename Base::StorageIndex StorageIndex;
446 enum { UpLo = _UpLo };
452 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
455 explicit PardisoLLT(
const MatrixType& matrix)
458 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
464 void getMatrix(
const MatrixType& matrix)
467 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
468 m_matrix.resize(matrix.rows(), matrix.cols());
469 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
470 m_matrix.makeCompressed();
492 template<
typename MatrixType,
int Options>
493 class PardisoLDLT :
public PardisoImpl< PardisoLDLT<MatrixType,Options> >
496 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
497 typedef typename Base::Scalar Scalar;
498 typedef typename Base::RealScalar RealScalar;
499 using Base::pardisoInit;
500 using Base::m_matrix;
501 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
505 typedef typename Base::StorageIndex StorageIndex;
512 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
515 explicit PardisoLDLT(
const MatrixType& matrix)
518 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
522 void getMatrix(
const MatrixType& matrix)
525 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
526 m_matrix.resize(matrix.rows(), matrix.cols());
527 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
528 m_matrix.makeCompressed();
534 #endif // EIGEN_PARDISOSUPPORT_H
Definition: Constants.h:222
Definition: Constants.h:206
const unsigned int RowMajorBit
Definition: Constants.h:61
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:38
Definition: Constants.h:434
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
Definition: Constants.h:439
Definition: Constants.h:432
Definition: Constants.h:204
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:430