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,
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, internal::convert_index<StorageIndex>(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 EIGEN_USING_STD_MATH(
abs);
196 bool symmetric =
abs(m_type) < 10;
207 m_iparm[10] = symmetric ? 0 : 1;
209 m_iparm[12] = symmetric ? 0 : 1;
221 m_iparm[27] = (
sizeof(RealScalar) == 4) ? 1 : 0;
226 memset(m_pt, 0,
sizeof(m_pt));
232 void manageErrorCode(
Index error)
const
248 mutable SparseMatrixType m_matrix;
250 bool m_analysisIsOk, m_factorizationIsOk;
251 StorageIndex m_type, m_msglvl;
252 mutable void *m_pt[64];
253 mutable ParameterType m_iparm;
254 mutable IntColVectorType m_perm;
259 template<
class Derived>
260 Derived& PardisoImpl<Derived>::compute(
const MatrixType& a)
263 eigen_assert(a.rows() == a.cols());
266 m_perm.setZero(m_size);
267 derived().getMatrix(a);
270 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size),
271 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
272 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
273 manageErrorCode(error);
274 m_analysisIsOk =
true;
275 m_factorizationIsOk =
true;
276 m_isInitialized =
true;
280 template<
class Derived>
281 Derived& PardisoImpl<Derived>::analyzePattern(
const MatrixType& a)
284 eigen_assert(m_size == a.cols());
287 m_perm.setZero(m_size);
288 derived().getMatrix(a);
291 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, internal::convert_index<StorageIndex>(m_size),
292 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
293 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
295 manageErrorCode(error);
296 m_analysisIsOk =
true;
297 m_factorizationIsOk =
false;
298 m_isInitialized =
true;
302 template<
class Derived>
303 Derived& PardisoImpl<Derived>::factorize(
const MatrixType& a)
305 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
306 eigen_assert(m_size == a.rows() && m_size == a.cols());
308 derived().getMatrix(a);
311 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, internal::convert_index<StorageIndex>(m_size),
312 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
313 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
315 manageErrorCode(error);
316 m_factorizationIsOk =
true;
320 template<
class Derived>
321 template<
typename BDerived,
typename XDerived>
322 void PardisoImpl<Derived>::_solve_impl(
const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x)
const
332 eigen_assert(m_size==b.rows());
335 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
347 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.derived().data());
348 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
351 if(rhs_ptr == x.derived().data())
354 rhs_ptr = tmp.data();
358 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, internal::convert_index<StorageIndex>(m_size),
359 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
360 m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl,
361 rhs_ptr, x.derived().data());
363 manageErrorCode(error);
384 template<
typename MatrixType>
385 class PardisoLU :
public PardisoImpl< PardisoLU<MatrixType> >
388 typedef PardisoImpl<PardisoLU> Base;
389 typedef typename Base::Scalar Scalar;
390 typedef typename Base::RealScalar RealScalar;
391 using Base::pardisoInit;
392 using Base::m_matrix;
393 friend class PardisoImpl< PardisoLU<MatrixType> >;
403 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
406 explicit PardisoLU(
const MatrixType& matrix)
409 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
413 void getMatrix(
const MatrixType& matrix)
416 m_matrix.makeCompressed();
439 template<
typename MatrixType,
int _UpLo>
440 class PardisoLLT :
public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
443 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
444 typedef typename Base::Scalar Scalar;
445 typedef typename Base::RealScalar RealScalar;
446 using Base::pardisoInit;
447 using Base::m_matrix;
448 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
452 typedef typename Base::StorageIndex StorageIndex;
453 enum { UpLo = _UpLo };
459 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
462 explicit PardisoLLT(
const MatrixType& matrix)
465 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
471 void getMatrix(
const MatrixType& matrix)
474 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
475 m_matrix.resize(matrix.rows(), matrix.cols());
476 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
477 m_matrix.makeCompressed();
502 template<
typename MatrixType,
int Options>
503 class PardisoLDLT :
public PardisoImpl< PardisoLDLT<MatrixType,Options> >
506 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
507 typedef typename Base::Scalar Scalar;
508 typedef typename Base::RealScalar RealScalar;
509 using Base::pardisoInit;
510 using Base::m_matrix;
511 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
515 typedef typename Base::StorageIndex StorageIndex;
522 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
525 explicit PardisoLDLT(
const MatrixType& matrix)
528 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
532 void getMatrix(
const MatrixType& matrix)
535 PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
536 m_matrix.resize(matrix.rows(), matrix.cols());
537 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
538 m_matrix.makeCompressed();
544 #endif // EIGEN_PARDISOSUPPORT_H