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>
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;
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>
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 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;
225 memset(m_pt, 0,
sizeof(m_pt));
231 void manageErrorCode(
Index error)
const 247 mutable SparseMatrixType m_matrix;
249 bool m_analysisIsOk, m_factorizationIsOk;
250 StorageIndex m_type, m_msglvl;
251 mutable void *m_pt[64];
252 mutable ParameterType m_iparm;
253 mutable IntColVectorType m_perm;
258 template<
class Derived>
259 Derived& PardisoImpl<Derived>::compute(
const MatrixType& a)
262 eigen_assert(a.rows() == a.cols());
265 m_perm.setZero(m_size);
266 derived().getMatrix(a);
269 error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size),
270 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
271 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, internal::convert_index<StorageIndex>(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, internal::convert_index<StorageIndex>(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>
331 eigen_assert(m_size==b.
rows());
346 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.
derived().data());
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, internal::convert_index<StorageIndex>(m_size),
358 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
359 m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl,
362 manageErrorCode(error);
383 template<
typename MatrixType>
384 class PardisoLU :
public PardisoImpl< PardisoLU<MatrixType> >
387 typedef PardisoImpl<PardisoLU> Base;
388 typedef typename Base::Scalar Scalar;
389 typedef typename Base::RealScalar RealScalar;
390 using Base::pardisoInit;
391 using Base::m_matrix;
392 friend class PardisoImpl<
PardisoLU<MatrixType> >;
402 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
405 explicit PardisoLU(
const MatrixType& matrix)
408 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
412 void getMatrix(
const MatrixType& matrix)
415 m_matrix.makeCompressed();
438 template<
typename MatrixType,
int _UpLo>
439 class PardisoLLT :
public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
442 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
443 typedef typename Base::Scalar Scalar;
444 typedef typename Base::RealScalar RealScalar;
445 using Base::pardisoInit;
446 using Base::m_matrix;
447 friend class PardisoImpl<
PardisoLLT<MatrixType,_UpLo> >;
451 typedef typename Base::StorageIndex StorageIndex;
452 enum { UpLo = _UpLo };
458 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
464 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
470 void getMatrix(
const MatrixType& matrix)
474 m_matrix.
resize(matrix.rows(), matrix.cols());
475 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
476 m_matrix.makeCompressed();
501 template<
typename MatrixType,
int Options>
502 class PardisoLDLT :
public PardisoImpl< PardisoLDLT<MatrixType,Options> >
505 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
506 typedef typename Base::Scalar Scalar;
507 typedef typename Base::RealScalar RealScalar;
508 using Base::pardisoInit;
509 using Base::m_matrix;
510 friend class PardisoImpl<
PardisoLDLT<MatrixType,Options> >;
514 typedef typename Base::StorageIndex StorageIndex;
521 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
524 explicit PardisoLDLT(
const MatrixType& matrix)
527 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
531 void getMatrix(
const MatrixType& matrix)
535 m_matrix.
resize(matrix.rows(), matrix.cols());
536 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
537 m_matrix.makeCompressed();
543 #endif // EIGEN_PARDISOSUPPORT_H
A base class for sparse solvers.
Definition: SparseSolverBase.h:67
Namespace containing all symbols from the Eigen library.
Definition: Core:287
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Derived & derived()
Definition: EigenBase.h:45
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:204
Index rows() const
Definition: EigenBase.h:59
Definition: Constants.h:222
Definition: Constants.h:434
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
Definition: Constants.h:439
Definition: Constants.h:206
Definition: Constants.h:432
const Scalar * data() const
Definition: PlainObjectBase.h:255
Definition: Eigen_Colamd.h:50
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:39
const int Dynamic
Definition: Constants.h:21
Index outerStride() const
Definition: DenseCoeffsBase.h:581
Index cols() const
Definition: EigenBase.h:62
void resize(Index newSize)
Definition: PermutationMatrix.h:136
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48