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 Index>
44 struct pardiso_run_selector
46 static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n,
void *a,
47 Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index 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 Index;
58 static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n,
void *a,
59 Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index 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::Index Index;
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::Index Index;
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::Index Index;
98 template<
class Derived>
101 typedef internal::pardiso_traits<Derived> Traits;
103 typedef typename Traits::MatrixType MatrixType;
104 typedef typename Traits::Scalar Scalar;
105 typedef typename Traits::RealScalar RealScalar;
106 typedef typename Traits::Index Index;
118 eigen_assert((
sizeof(Index) >=
sizeof(_INTEGER_t) &&
sizeof(Index) <= 8) &&
"Non-supported index type");
121 m_initialized =
false;
129 inline Index cols()
const {
return m_size; }
130 inline Index rows()
const {
return m_size; }
139 eigen_assert(m_initialized &&
"Decomposition is not initialized.");
146 ParameterType& pardisoParameterArray()
157 Derived& analyzePattern(
const MatrixType& matrix);
165 Derived& factorize(
const MatrixType& matrix);
167 Derived& compute(
const MatrixType& matrix);
173 template<
typename Rhs>
174 inline const internal::solve_retval<PardisoImpl, Rhs>
177 eigen_assert(m_initialized &&
"Pardiso solver is not initialized.");
178 eigen_assert(rows()==b.rows()
179 &&
"PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
180 return internal::solve_retval<PardisoImpl, Rhs>(*
this, b.derived());
187 template<
typename Rhs>
188 inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
191 eigen_assert(m_initialized &&
"Pardiso solver is not initialized.");
192 eigen_assert(rows()==b.
rows()
193 &&
"PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
194 return internal::sparse_solve_retval<PardisoImpl, Rhs>(*
this, b.
derived());
199 return *
static_cast<Derived*
>(
this);
201 const Derived& derived()
const 203 return *
static_cast<const Derived*
>(
this);
206 template<
typename BDerived,
typename XDerived>
210 void pardisoRelease()
214 internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
215 m_iparm.data(), m_msglvl, 0, 0);
219 void pardisoInit(
int type)
222 bool symmetric = std::abs(m_type) < 10;
233 m_iparm[10] = symmetric ? 0 : 1;
235 m_iparm[12] = symmetric ? 0 : 1;
247 m_iparm[27] = (
sizeof(RealScalar) == 4) ? 1 : 0;
252 memset(m_pt, 0,
sizeof(m_pt));
258 void manageErrorCode(Index error)
274 mutable SparseMatrixType m_matrix;
276 bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
277 Index m_type, m_msglvl;
278 mutable void *m_pt[64];
279 mutable ParameterType m_iparm;
280 mutable IntColVectorType m_perm;
284 PardisoImpl(PardisoImpl &) {}
287 template<
class Derived>
288 Derived& PardisoImpl<Derived>::compute(
const MatrixType& a)
291 eigen_assert(a.rows() == a.cols());
294 memset(m_pt, 0,
sizeof(m_pt));
295 m_perm.setZero(m_size);
296 derived().getMatrix(a);
299 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
300 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
301 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
303 manageErrorCode(error);
304 m_analysisIsOk =
true;
305 m_factorizationIsOk =
true;
306 m_initialized =
true;
310 template<
class Derived>
311 Derived& PardisoImpl<Derived>::analyzePattern(
const MatrixType& a)
314 eigen_assert(m_size == a.cols());
317 memset(m_pt, 0,
sizeof(m_pt));
318 m_perm.setZero(m_size);
319 derived().getMatrix(a);
322 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
323 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
324 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
326 manageErrorCode(error);
327 m_analysisIsOk =
true;
328 m_factorizationIsOk =
false;
329 m_initialized =
true;
333 template<
class Derived>
334 Derived& PardisoImpl<Derived>::factorize(
const MatrixType& a)
336 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
337 eigen_assert(m_size == a.rows() && m_size == a.cols());
339 derived().getMatrix(a);
342 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
343 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
344 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
346 manageErrorCode(error);
347 m_factorizationIsOk =
true;
352 template<
typename BDerived,
typename XDerived>
359 Index nrhs = Index(b.cols());
360 eigen_assert(m_size==b.rows());
363 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
375 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.derived().data());
379 if(rhs_ptr == x.derived().data())
382 rhs_ptr = tmp.
data();
386 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
387 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
388 m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
389 rhs_ptr, x.derived().data());
409 template<
typename MatrixType>
410 class PardisoLU :
public PardisoImpl< PardisoLU<MatrixType> >
413 typedef PardisoImpl< PardisoLU<MatrixType> > Base;
414 typedef typename Base::Scalar Scalar;
415 typedef typename Base::RealScalar RealScalar;
416 using Base::pardisoInit;
417 using Base::m_matrix;
418 friend class PardisoImpl<
PardisoLU<MatrixType> >;
428 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
434 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
438 void getMatrix(
const MatrixType& matrix)
464 template<
typename MatrixType,
int _UpLo>
465 class PardisoLLT :
public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
468 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
469 typedef typename Base::Scalar Scalar;
470 typedef typename Base::Index Index;
471 typedef typename Base::RealScalar RealScalar;
472 using Base::pardisoInit;
473 using Base::m_matrix;
474 friend class PardisoImpl<
PardisoLLT<MatrixType,_UpLo> >;
478 enum { UpLo = _UpLo };
485 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
491 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
497 void getMatrix(
const MatrixType& matrix)
501 m_matrix.
resize(matrix.rows(), matrix.cols());
502 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
528 template<
typename MatrixType,
int Options>
529 class PardisoLDLT :
public PardisoImpl< PardisoLDLT<MatrixType,Options> >
532 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
533 typedef typename Base::Scalar Scalar;
534 typedef typename Base::Index Index;
535 typedef typename Base::RealScalar RealScalar;
536 using Base::pardisoInit;
537 using Base::m_matrix;
538 friend class PardisoImpl<
PardisoLDLT<MatrixType,Options> >;
549 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
552 PardisoLDLT(
const MatrixType& matrix)
555 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
559 void getMatrix(
const MatrixType& matrix)
563 m_matrix.
resize(matrix.rows(), matrix.cols());
564 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
568 PardisoLDLT(PardisoLDLT& ) {}
573 template<
typename _Derived,
typename Rhs>
574 struct solve_retval<PardisoImpl<_Derived>, Rhs>
575 : solve_retval_base<PardisoImpl<_Derived>, Rhs>
577 typedef PardisoImpl<_Derived> Dec;
578 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
580 template<typename Dest>
void evalTo(Dest& dst)
const 582 dec()._solve(rhs(),dst);
586 template<
typename Derived,
typename Rhs>
587 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
588 : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
590 typedef PardisoImpl<Derived> Dec;
591 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
593 template<typename Dest>
void evalTo(Dest& dst)
const 595 this->defaultEvalTo(dst);
603 #endif // EIGEN_PARDISOSUPPORT_H const Scalar * data() const
Definition: PlainObjectBase.h:212
Definition: Constants.h:378
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Definition: Constants.h:185
Definition: Constants.h:167
void resize(Index newSize)
Definition: PermutationMatrix.h:142
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:38
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:239
Derived & derived()
Definition: EigenBase.h:34
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
Definition: Constants.h:383
Definition: Eigen_Colamd.h:50
Definition: Constants.h:169
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:39
Definition: Constants.h:376
const unsigned int RowMajorBit
Definition: Constants.h:53
Index rows() const
Definition: SparseMatrixBase.h:160
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48