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>
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);
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);
69 template<
typename _MatrixType>
73 typedef typename _MatrixType::Scalar
Scalar;
75 typedef typename _MatrixType::Index
Index;
78 template<
typename _MatrixType,
int Options>
82 typedef typename _MatrixType::Scalar
Scalar;
84 typedef typename _MatrixType::Index
Index;
87 template<
typename _MatrixType,
int Options>
91 typedef typename _MatrixType::Scalar
Scalar;
93 typedef typename _MatrixType::Index
Index;
98 template<
class Derived>
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.");
157 Derived& analyzePattern(
const MatrixType& matrix);
165 Derived& factorize(
const MatrixType& matrix);
167 Derived& compute(
const MatrixType& matrix);
173 template<
typename Rhs>
177 eigen_assert(m_initialized &&
"Pardiso solver is not initialized.");
179 &&
"PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
187 template<
typename Rhs>
191 eigen_assert(m_initialized &&
"Pardiso solver is not initialized.");
193 &&
"PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
199 return *
static_cast<Derived*
>(
this);
203 return *
static_cast<const Derived*
>(
this);
206 template<
typename BDerived,
typename XDerived>
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);
222 bool symmetric =
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;
275 mutable void *m_pt[64];
284 template<
class Derived>
291 memset(m_pt, 0,
sizeof(m_pt));
292 m_perm.setZero(m_size);
293 derived().getMatrix(a);
297 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
298 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
300 manageErrorCode(error);
301 m_analysisIsOk =
true;
302 m_factorizationIsOk =
true;
303 m_initialized =
true;
307 template<
class Derived>
314 memset(m_pt, 0,
sizeof(m_pt));
315 m_perm.setZero(m_size);
316 derived().getMatrix(a);
320 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
321 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
323 manageErrorCode(error);
324 m_analysisIsOk =
true;
325 m_factorizationIsOk =
false;
326 m_initialized =
true;
330 template<
class Derived>
333 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
336 derived().getMatrix(a);
340 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
341 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
343 manageErrorCode(error);
344 m_factorizationIsOk =
true;
349 template<
typename BDerived,
typename XDerived>
360 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
372 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.derived().data());
376 if(rhs_ptr == x.derived().data())
379 rhs_ptr = tmp.
data();
384 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
385 m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
386 rhs_ptr, x.derived().data());
404 template<
typename MatrixType>
411 using Base::pardisoInit;
412 using Base::m_matrix;
423 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
429 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
456 template<
typename MatrixType,
int _UpLo>
464 using Base::pardisoInit;
465 using Base::m_matrix;
470 enum { UpLo = _UpLo };
477 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
483 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
493 m_matrix.
resize(matrix.rows(), matrix.cols());
494 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
517 template<
typename MatrixType,
int Options>
525 using Base::pardisoInit;
526 using Base::m_matrix;
538 pardisoInit(Base::ScalarIsComplex ? (
bool(
Options&
Symmetric) ? 6 : -4 ) : -2);
544 pardisoInit(Base::ScalarIsComplex ? (
bool(
Options&
Symmetric) ? 6 : -4 ) : -2);
552 m_matrix.
resize(matrix.rows(), matrix.cols());
553 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
562 template<
typename _Derived,
typename Rhs>
569 template<typename Dest>
void evalTo(Dest& dst)
const 571 dec()._solve(
rhs(),dst);
575 template<
typename Derived,
typename Rhs>
582 template<typename Dest>
void evalTo(Dest& dst)
const 584 this->defaultEvalTo(dst);
592 #endif // EIGEN_PARDISOSUPPORT_H
static Index run(_MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a, Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
_MatrixType::Scalar Scalar
PardisoImpl< PardisoLU< MatrixType > > Base
PardisoImpl< _Derived > Dec
Derived & compute(const MatrixType &matrix)
Matrix< Index, MatrixType::RowsAtCompileTime, 1 > IntColVectorType
Derived & factorize(const MatrixType &matrix)
PardisoLDLT(const MatrixType &matrix)
PardisoLDLT(PardisoLDLT &)
iterative scaling algorithm to equilibrate rows and column norms in matrices
const Derived & derived() const
PardisoImpl< PardisoLLT< MatrixType, _UpLo > > Base
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
const internal::solve_retval< PardisoImpl, Rhs > solve(const MatrixBase< Rhs > &b) const
void getMatrix(const MatrixType &matrix)
_MatrixType::RealScalar RealScalar
const unsigned int RowMajorBit
void manageErrorCode(Index error)
Traits::RealScalar RealScalar
PardisoImpl< Derived > Dec
PardisoImpl(PardisoImpl &)
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
Base::RealScalar RealScalar
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
SparseMatrix< Scalar, RowMajor, Index > SparseMatrixType
void pardisoInit(int type)
Base class of any sparse matrices or sparse expressions.
Derived & analyzePattern(const MatrixType &matrix)
_MatrixType::RealScalar RealScalar
const internal::sparse_solve_retval< PardisoImpl, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
Provides a generic way to set and pass user-specified options.
Base::RealScalar RealScalar
static Index run(_MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a, Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
A sparse direct LU factorization and solver based on the PARDISO library.
Matrix< Scalar, Dynamic, 1 > VectorType
#define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType, Rhs)
void rhs(const real_t *x, real_t *f)
EIGEN_STRONG_INLINE const Scalar * data() const
internal::pardiso_traits< Derived > Traits
const Derived & derived() const
ComputationInfo info() const
Reports whether previous computation was successful.
void getMatrix(const MatrixType &matrix)
Matrix< Index, 1, MatrixType::ColsAtCompileTime > IntRowVectorType
Base::RealScalar RealScalar
PardisoImpl< PardisoLDLT< MatrixType, Options > > Base
bool _solve(const MatrixBase< BDerived > &b, MatrixBase< XDerived > &x) const
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
ParameterType & pardisoParameterArray()
void getMatrix(const MatrixType &matrix)
Array< Index, 64, 1, DontAlign > ParameterType
Traits::MatrixType MatrixType
void resize(Index newSize)
Base class for all dense matrices, vectors, and expressions.
SparseMatrixType m_matrix
_MatrixType::Scalar Scalar
_MatrixType::RealScalar RealScalar
PardisoLLT(const MatrixType &matrix)
PardisoLU(const MatrixType &matrix)
_MatrixType::Scalar Scalar