32 #ifndef EIGEN_PARDISOSUPPORT_H
33 #define EIGEN_PARDISOSUPPORT_H
37 template<
typename _MatrixType>
class PardisoLU;
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;
106 typedef typename Traits::Index
Index;
172 template<
typename Rhs>
173 inline const internal::solve_retval<PardisoImpl, Rhs>
178 &&
"PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
179 return internal::solve_retval<PardisoImpl, Rhs>(*
this, b.derived());
186 template<
typename Rhs>
187 inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
192 &&
"PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
193 return internal::sparse_solve_retval<PardisoImpl, Rhs>(*
this, b.
derived());
198 return *
static_cast<Derived*
>(
this);
202 return *
static_cast<const Derived*
>(
this);
205 template<
typename BDerived,
typename XDerived>
209 template<
typename Rhs,
typename DestScalar,
int DestOptions,
typename DestIndex>
215 static const int NbColsAtOnce = 4;
216 int rhsCols = b.cols();
222 for(
int k=0; k<rhsCols; k+=NbColsAtOnce)
224 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
225 tmp_rhs.
leftCols(actualCols) = b.middleCols(k,actualCols);
236 internal::pardiso_run_selector<Index>::run(
m_pt, 1, 1,
m_type, -1,
m_size, 0, 0, 0,
m_perm.
data(), 0,
255 m_iparm[10] = symmetric ? 0 : 1;
257 m_iparm[12] = symmetric ? 0 : 1;
306 template<
class Derived>
313 memset(m_pt, 0,
sizeof(m_pt));
314 m_perm.setZero(m_size);
315 derived().getMatrix(a);
318 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
319 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
320 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
322 manageErrorCode(error);
323 m_analysisIsOk =
true;
324 m_factorizationIsOk =
true;
325 m_initialized =
true;
329 template<
class Derived>
336 memset(m_pt, 0,
sizeof(m_pt));
337 m_perm.setZero(m_size);
338 derived().getMatrix(a);
341 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
342 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
343 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
345 manageErrorCode(error);
346 m_analysisIsOk =
true;
347 m_factorizationIsOk =
false;
348 m_initialized =
true;
352 template<
class Derived>
355 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
358 derived().getMatrix(a);
361 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
362 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
363 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
365 manageErrorCode(error);
366 m_factorizationIsOk =
true;
371 template<
typename BDerived,
typename XDerived>
382 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
394 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.derived().data());
398 if(rhs_ptr == x.derived().data())
401 rhs_ptr = tmp.
data();
405 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
406 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
407 m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
408 rhs_ptr, x.derived().data());
426 template<
typename MatrixType>
478 template<
typename MatrixType,
int _UpLo>
516 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().
twistedBy(p_null);
539 template<
typename MatrixType,
int Options>
575 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().
twistedBy(p_null);
584 template<
typename _Derived,
typename Rhs>
586 : solve_retval_base<PardisoImpl<_Derived>, Rhs>
588 typedef PardisoImpl<_Derived> Dec;
591 template<typename Dest>
void evalTo(Dest& dst)
const
593 dec()._solve(rhs(),dst);
597 template<
typename Derived,
typename Rhs>
598 struct sparse_solve_retval<
PardisoImpl<Derived>, Rhs>
599 : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
601 typedef PardisoImpl<Derived> Dec;
604 template<typename Dest>
void evalTo(Dest& dst)
const
606 dec().derived()._solve_sparse(rhs(),dst);
614 #endif // EIGEN_PARDISOSUPPORT_H