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>
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;
103 typedef typename Traits::MatrixType MatrixType;
104 typedef typename Traits::Scalar Scalar;
105 typedef typename Traits::RealScalar RealScalar;
106 typedef typename Traits::Index Index;
107 typedef SparseMatrix<Scalar,RowMajor,Index> SparseMatrixType;
108 typedef Matrix<Scalar,Dynamic,1> VectorType;
109 typedef Matrix<Index, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
110 typedef Matrix<Index, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
111 typedef Array<Index,64,1,DontAlign> ParameterType;
113 ScalarIsComplex = NumTraits<Scalar>::IsComplex
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.");
146 ParameterType& pardisoParameterArray()
157 Derived& analyzePattern(
const MatrixType& matrix);
165 Derived& factorize(
const MatrixType& matrix);
167 Derived& compute(
const MatrixType& matrix);
173 template<
typename Rhs>
174 inline const internal::solve_retval<PardisoImpl, Rhs>
175 solve(
const MatrixBase<Rhs>& b)
const
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");
180 return internal::solve_retval<PardisoImpl, Rhs>(*
this, b.derived());
187 template<
typename Rhs>
188 inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
189 solve(
const SparseMatrixBase<Rhs>& b)
const
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");
194 return internal::sparse_solve_retval<PardisoImpl, Rhs>(*
this, b.derived());
199 return *
static_cast<Derived*
>(
this);
201 const Derived& derived()
const
203 return *
static_cast<const Derived*
>(
this);
206 template<
typename BDerived,
typename XDerived>
207 bool _solve(
const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x)
const;
210 template<
typename Rhs,
typename DestScalar,
int DestOptions,
typename DestIndex>
211 void _solve_sparse(
const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest)
const
213 eigen_assert(m_size==b.rows());
216 static const int NbColsAtOnce = 4;
217 int rhsCols = b.cols();
223 for(
int k=0; k<rhsCols; k+=NbColsAtOnce)
225 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
226 tmp_rhs.leftCols(actualCols) = b.middleCols(k,actualCols);
227 tmp_res.leftCols(actualCols) = derived().solve(tmp_rhs.leftCols(actualCols));
228 dest.middleCols(k,actualCols) = tmp_res.leftCols(actualCols).sparseView();
233 void pardisoRelease()
237 internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
238 m_iparm.data(), m_msglvl, 0, 0);
242 void pardisoInit(
int type)
245 bool symmetric = abs(m_type) < 10;
256 m_iparm[10] = symmetric ? 0 : 1;
258 m_iparm[12] = symmetric ? 0 : 1;
270 m_iparm[27] = (
sizeof(RealScalar) == 4) ? 1 : 0;
278 void manageErrorCode(Index error)
294 mutable SparseMatrixType m_matrix;
296 bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
297 Index m_type, m_msglvl;
298 mutable void *m_pt[64];
299 mutable ParameterType m_iparm;
300 mutable IntColVectorType m_perm;
304 PardisoImpl(PardisoImpl &) {}
307 template<
class Derived>
308 Derived& PardisoImpl<Derived>::compute(
const MatrixType& a)
311 eigen_assert(a.rows() == a.cols());
314 memset(m_pt, 0,
sizeof(m_pt));
315 m_perm.setZero(m_size);
316 derived().getMatrix(a);
319 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
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 =
true;
326 m_initialized =
true;
330 template<
class Derived>
331 Derived& PardisoImpl<Derived>::analyzePattern(
const MatrixType& a)
334 eigen_assert(m_size == a.cols());
337 memset(m_pt, 0,
sizeof(m_pt));
338 m_perm.setZero(m_size);
339 derived().getMatrix(a);
342 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
343 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
344 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
346 manageErrorCode(error);
347 m_analysisIsOk =
true;
348 m_factorizationIsOk =
false;
349 m_initialized =
true;
353 template<
class Derived>
354 Derived& PardisoImpl<Derived>::factorize(
const MatrixType& a)
356 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
357 eigen_assert(m_size == a.rows() && m_size == a.cols());
359 derived().getMatrix(a);
362 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
363 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
364 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
366 manageErrorCode(error);
367 m_factorizationIsOk =
true;
372 template<
typename BDerived,
typename XDerived>
373 bool PardisoImpl<Base>::_solve(
const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x)
const
379 Index nrhs = Index(b.cols());
380 eigen_assert(m_size==b.rows());
381 eigen_assert(((MatrixBase<BDerived>::Flags &
RowMajorBit) == 0 || nrhs == 1) &&
"Row-major right hand sides are not supported");
382 eigen_assert(((MatrixBase<XDerived>::Flags &
RowMajorBit) == 0 || nrhs == 1) &&
"Row-major matrices of unknowns are not supported");
383 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
395 Scalar* rhs_ptr =
const_cast<Scalar*
>(b.derived().data());
396 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
399 if(rhs_ptr == x.derived().data())
402 rhs_ptr = tmp.data();
406 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
407 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
408 m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
409 rhs_ptr, x.derived().data());
427 template<
typename MatrixType>
428 class PardisoLU :
public PardisoImpl< PardisoLU<MatrixType> >
431 typedef PardisoImpl< PardisoLU<MatrixType> > Base;
432 typedef typename Base::Scalar Scalar;
433 typedef typename Base::RealScalar RealScalar;
434 using Base::pardisoInit;
435 using Base::m_matrix;
436 friend class PardisoImpl< PardisoLU<MatrixType> >;
446 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
449 PardisoLU(
const MatrixType& matrix)
452 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
456 void getMatrix(
const MatrixType& matrix)
462 PardisoLU(PardisoLU& ) {}
479 template<
typename MatrixType,
int _UpLo>
480 class PardisoLLT :
public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
483 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
484 typedef typename Base::Scalar Scalar;
485 typedef typename Base::Index Index;
486 typedef typename Base::RealScalar RealScalar;
487 using Base::pardisoInit;
488 using Base::m_matrix;
489 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
493 enum { UpLo = _UpLo };
500 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
503 PardisoLLT(
const MatrixType& matrix)
506 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
512 void getMatrix(
const MatrixType& matrix)
515 PermutationMatrix<Dynamic,Dynamic,Index> p_null;
516 m_matrix.resize(matrix.rows(), matrix.cols());
517 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
521 PardisoLLT(PardisoLLT& ) {}
540 template<
typename MatrixType,
int Options>
541 class PardisoLDLT :
public PardisoImpl< PardisoLDLT<MatrixType,Options> >
544 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
545 typedef typename Base::Scalar Scalar;
546 typedef typename Base::Index Index;
547 typedef typename Base::RealScalar RealScalar;
548 using Base::pardisoInit;
549 using Base::m_matrix;
550 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
561 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
564 PardisoLDLT(
const MatrixType& matrix)
567 pardisoInit(Base::ScalarIsComplex ? (
bool(Options&
Symmetric) ? 6 : -4 ) : -2);
571 void getMatrix(
const MatrixType& matrix)
574 PermutationMatrix<Dynamic,Dynamic,Index> p_null;
575 m_matrix.resize(matrix.rows(), matrix.cols());
576 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
580 PardisoLDLT(PardisoLDLT& ) {}
585 template<
typename _Derived,
typename Rhs>
586 struct solve_retval<PardisoImpl<_Derived>, Rhs>
587 : solve_retval_base<PardisoImpl<_Derived>, Rhs>
589 typedef PardisoImpl<_Derived> Dec;
590 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
592 template<typename Dest>
void evalTo(Dest& dst)
const
594 dec()._solve(rhs(),dst);
598 template<
typename Derived,
typename Rhs>
599 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
600 : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
602 typedef PardisoImpl<Derived> Dec;
603 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
605 template<typename Dest>
void evalTo(Dest& dst)
const
607 dec().derived()._solve_sparse(rhs(),dst);
615 #endif // EIGEN_PARDISOSUPPORT_H