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;
126 : m_analysisIsOk(false), m_factorizationIsOk(false)
128 eigen_assert((
sizeof(StorageIndex) >=
sizeof(_INTEGER_t) &&
sizeof(StorageIndex) <= 8) &&
"Non-supported index type");
131 m_isInitialized =
false;
139 inline Index cols()
const {
return m_size; }
140 inline Index rows()
const {
return m_size; }
149 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
156 ParameterType& pardisoParameterArray()
167 Derived& analyzePattern(
const MatrixType& matrix);
175 Derived& factorize(
const MatrixType& matrix);
177 Derived& compute(
const MatrixType& matrix);
179 template<
typename Rhs,
typename Dest>
183 void pardisoRelease()
187 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,
188 m_iparm.data(), m_msglvl, NULL, NULL);
189 m_isInitialized =
false;
193 void pardisoInit(
int type)
196 bool symmetric = std::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>
332 eigen_assert(m_size==b.
rows());
347 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.
derived().data());
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,
363 manageErrorCode(error);
384 template<
typename MatrixType>
385 class PardisoLU :
public PardisoImpl< PardisoLU<MatrixType> >
388 typedef PardisoImpl<PardisoLU> Base;
389 using Base::pardisoInit;
390 using Base::m_matrix;
391 friend class PardisoImpl<
PardisoLU<MatrixType> >;
395 typedef typename Base::Scalar Scalar;
396 typedef typename Base::RealScalar RealScalar;
404 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
407 explicit PardisoLU(
const MatrixType& matrix)
410 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
414 void getMatrix(
const MatrixType& matrix)
417 m_matrix.makeCompressed();
440 template<
typename MatrixType,
int _UpLo>
441 class PardisoLLT :
public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
444 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
445 using Base::pardisoInit;
446 using Base::m_matrix;
447 friend class PardisoImpl<
PardisoLLT<MatrixType,_UpLo> >;
451 typedef typename Base::Scalar Scalar;
452 typedef typename Base::RealScalar RealScalar;
453 typedef typename Base::StorageIndex StorageIndex;
454 enum { UpLo = _UpLo };
460 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
466 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
472 void getMatrix(
const MatrixType& matrix)
476 m_matrix.
resize(matrix.rows(), matrix.cols());
477 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
478 m_matrix.makeCompressed();
503 template<
typename MatrixType,
int Options>
504 class PardisoLDLT :
public PardisoImpl< PardisoLDLT<MatrixType,Options> >
507 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
508 using Base::pardisoInit;
509 using Base::m_matrix;
510 friend class PardisoImpl<
PardisoLDLT<MatrixType,Options> >;
514 typedef typename Base::Scalar Scalar;
515 typedef typename Base::RealScalar RealScalar;
516 typedef typename Base::StorageIndex StorageIndex;
523 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
526 explicit PardisoLDLT(
const MatrixType& matrix)
529 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
533 void getMatrix(
const MatrixType& matrix)
537 m_matrix.
resize(matrix.rows(), matrix.cols());
538 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
539 m_matrix.makeCompressed();
545 #endif // EIGEN_PARDISOSUPPORT_H EIGEN_CONSTEXPR Index outerStride() const EIGEN_NOEXCEPT
Definition: DenseCoeffsBase.h:585
A base class for sparse solvers.
Definition: SparseSolverBase.h:67
Namespace containing all symbols from the Eigen library.
Definition: Core:141
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Derived & derived()
Definition: EigenBase.h:46
const unsigned int RowMajorBit
Definition: Constants.h:66
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:38
Definition: Constants.h:209
Definition: Constants.h:227
Definition: Constants.h:444
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
Definition: Constants.h:449
Definition: Constants.h:211
Definition: Constants.h:442
const Scalar * data() const
Definition: PlainObjectBase.h:247
Definition: Eigen_Colamd.h:50
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:39
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
const int Dynamic
Definition: Constants.h:22
void resize(Index newSize)
Definition: PermutationMatrix.h:125
ComputationInfo
Definition: Constants.h:440
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48