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);
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);
69 template<
typename _MatrixType>
72 typedef _MatrixType MatrixType;
74 typedef typename _MatrixType::RealScalar RealScalar;
75 typedef typename _MatrixType::Index Index;
78 template<
typename _MatrixType,
int Options>
81 typedef _MatrixType MatrixType;
83 typedef typename _MatrixType::RealScalar RealScalar;
84 typedef typename _MatrixType::Index Index;
87 template<
typename _MatrixType,
int Options>
90 typedef _MatrixType MatrixType;
92 typedef typename _MatrixType::RealScalar RealScalar;
93 typedef typename _MatrixType::Index Index;
98 template<
class Derived>
103 typedef typename Traits::MatrixType MatrixType;
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.");
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.");
178 eigen_assert(rows()==b.rows()
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.");
192 eigen_assert(rows()==b.
rows()
193 &&
"PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
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;
255 void manageErrorCode(Index error)
271 mutable SparseMatrixType m_matrix;
273 bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
274 Index m_type, m_msglvl;
275 mutable void *m_pt[64];
276 mutable ParameterType m_iparm;
277 mutable IntColVectorType m_perm;
284 template<
class Derived>
288 eigen_assert(a.rows() == a.cols());
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>
311 eigen_assert(m_size == a.cols());
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()");
334 eigen_assert(m_size == a.rows() && m_size == a.cols());
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>
356 Index nrhs = Index(b.cols());
357 eigen_assert(m_size==b.rows());
359 eigen_assert(((
MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) &&
"Row-major matrices of unknowns are not supported");
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>
410 typedef typename Base::RealScalar RealScalar;
411 using Base::pardisoInit;
412 using Base::m_matrix;
423 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
426 PardisoLU(
const MatrixType& matrix)
429 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
433 void getMatrix(
const MatrixType& matrix)
439 PardisoLU(PardisoLU& ) {}
456 template<
typename MatrixType,
int _UpLo>
462 typedef typename Base::Index Index;
463 typedef typename Base::RealScalar RealScalar;
464 using Base::pardisoInit;
465 using Base::m_matrix;
466 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
470 enum { UpLo = _UpLo };
477 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
480 PardisoLLT(
const MatrixType& matrix)
483 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
489 void getMatrix(
const MatrixType& matrix)
493 m_matrix.
resize(matrix.rows(), matrix.cols());
494 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
498 PardisoLLT(PardisoLLT& ) {}
517 template<
typename MatrixType,
int Options>
523 typedef typename Base::Index Index;
524 typedef typename Base::RealScalar RealScalar;
525 using Base::pardisoInit;
526 using Base::m_matrix;
527 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
538 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
541 PardisoLDLT(
const MatrixType& matrix)
544 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
548 void getMatrix(
const MatrixType& matrix)
552 m_matrix.
resize(matrix.rows(), matrix.cols());
553 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
557 PardisoLDLT(PardisoLDLT& ) {}
562 template<
typename _Derived,
typename Rhs>
567 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
569 template<
typename Dest>
void evalTo(Dest& dst)
const 571 dec()._solve(rhs(),dst);
575 template<
typename Derived,
typename Rhs>
580 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
582 template<
typename Dest>
void evalTo(Dest& dst)
const 584 this->defaultEvalTo(dst);
592 #endif // EIGEN_PARDISOSUPPORT_H Definition: gtest_unittest.cc:5031
Used to support symmetric, non-selfadjoint, complex matrices.
Definition: Constants.h:185
Definition: ForwardDeclarations.h:124
Derived & factorize(const MatrixType &matrix)
Performs a numeric decomposition of matrix.
Definition: PardisoSupport.h:331
EIGEN_STRONG_INLINE const Scalar * data() const
Definition: PlainObjectBase.h:212
const internal::solve_retval< PardisoImpl, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PardisoSupport.h:175
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: TestIMU_Common.h:87
Definition: SparseSolve.h:18
Holds information about the various numeric (i.e.
Definition: NumTraits.h:88
const unsigned int RowMajorBit
for a matrix, this means that the storage order is row-major.
Definition: Constants.h:53
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:38
View matrix as an upper triangular matrix.
Definition: Constants.h:169
Index rows() const
Definition: SparseMatrixBase.h:159
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:239
Derived & analyzePattern(const MatrixType &matrix)
Performs a symbolic decomposition on the sparcity of matrix.
Definition: PardisoSupport.h:308
The provided data did not satisfy the prerequisites.
Definition: Constants.h:378
const internal::sparse_solve_retval< PardisoImpl, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
Definition: PardisoSupport.h:189
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
The inputs are invalid, or the algorithm has been improperly called.
Definition: Constants.h:383
Definition: PardisoSupport.h:44
Computation was successful.
Definition: Constants.h:376
Definition: PardisoSupport.h:99
Definition: BandTriangularSolver.h:13
Definition: SparseSolve.h:17
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: PardisoSupport.h:137
Definition: ForwardDeclarations.h:125
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:39
View matrix as a lower triangular matrix.
Definition: Constants.h:167
ParameterType & pardisoParameterArray()
Definition: PardisoSupport.h:146
void resize(Index newSize)
Resizes to given size.
Definition: PermutationMatrix.h:142
ComputationInfo
Enum for reporting the status of a computation.
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Definition: osvr_print_tree.cpp:52
double Scalar
Common scalar type.
Definition: FlexibleKalmanBase.h:48
Definition: PardisoSupport.h:67