25 #ifndef EIGEN_PASTIXSUPPORT_H
26 #define EIGEN_PASTIXSUPPORT_H
38 template<
typename _MatrixType,
bool IsStrSym = false>
class PastixLU;
39 template<
typename _MatrixType,
int Options>
class PastixLLT;
40 template<
typename _MatrixType,
int Options>
class PastixLDLT;
45 template<
class Pastix>
struct pastix_traits;
47 template<
typename _MatrixType>
48 struct pastix_traits< PastixLU<_MatrixType> >
50 typedef _MatrixType MatrixType;
51 typedef typename _MatrixType::Scalar Scalar;
52 typedef typename _MatrixType::RealScalar RealScalar;
53 typedef typename _MatrixType::Index Index;
56 template<
typename _MatrixType,
int Options>
57 struct pastix_traits< PastixLLT<_MatrixType,Options> >
59 typedef _MatrixType MatrixType;
60 typedef typename _MatrixType::Scalar Scalar;
61 typedef typename _MatrixType::RealScalar RealScalar;
62 typedef typename _MatrixType::Index Index;
65 template<
typename _MatrixType,
int Options>
66 struct pastix_traits< PastixLDLT<_MatrixType,Options> >
68 typedef _MatrixType MatrixType;
69 typedef typename _MatrixType::Scalar Scalar;
70 typedef typename _MatrixType::RealScalar RealScalar;
71 typedef typename _MatrixType::Index Index;
74 void eigen_pastix(pastix_data_t **pastix_data,
int pastix_comm,
int n,
int *ptr,
int *idx,
float *vals,
int *perm,
int * invp,
float *x,
int nbrhs,
int *iparm,
double *dparm)
76 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
77 if (nbrhs == 0) {x = NULL; nbrhs=1;}
78 s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
81 void eigen_pastix(pastix_data_t **pastix_data,
int pastix_comm,
int n,
int *ptr,
int *idx,
double *vals,
int *perm,
int * invp,
double *x,
int nbrhs,
int *iparm,
double *dparm)
83 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
84 if (nbrhs == 0) {x = NULL; nbrhs=1;}
85 d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
88 void eigen_pastix(pastix_data_t **pastix_data,
int pastix_comm,
int n,
int *ptr,
int *idx, std::complex<float> *vals,
int *perm,
int * invp, std::complex<float> *x,
int nbrhs,
int *iparm,
double *dparm)
90 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
91 if (nbrhs == 0) {x = NULL; nbrhs=1;}
92 c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm);
95 void eigen_pastix(pastix_data_t **pastix_data,
int pastix_comm,
int n,
int *ptr,
int *idx, std::complex<double> *vals,
int *perm,
int * invp, std::complex<double> *x,
int nbrhs,
int *iparm,
double *dparm)
97 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
98 if (nbrhs == 0) {x = NULL; nbrhs=1;}
99 z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm);
103 template <
typename MatrixType>
106 if ( !(mat.outerIndexPtr()[0]) )
109 for(i = 0; i <= mat.rows(); ++i)
110 ++mat.outerIndexPtr()[i];
111 for(i = 0; i < mat.nonZeros(); ++i)
112 ++mat.innerIndexPtr()[i];
117 template <
typename MatrixType>
121 if ( mat.outerIndexPtr()[0] == 1 )
124 for(i = 0; i <= mat.rows(); ++i)
125 --mat.outerIndexPtr()[i];
126 for(i = 0; i < mat.nonZeros(); ++i)
127 --mat.innerIndexPtr()[i];
134 template <
class Derived>
138 typedef typename internal::pastix_traits<Derived>::MatrixType
_MatrixType;
140 typedef typename MatrixType::Scalar
Scalar;
142 typedef typename MatrixType::Index
Index;
162 template<
typename Rhs>
163 inline const internal::solve_retval<PastixBase, Rhs>
168 &&
"PastixBase::solve(): invalid number of rows of the right hand side matrix b");
169 return internal::solve_retval<PastixBase, Rhs>(*
this, b.derived());
172 template<
typename Rhs,
typename Dest>
176 template<
typename Rhs,
typename DestScalar,
int DestOptions,
typename DestIndex>
179 eigen_assert(
m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
183 static const int NbColsAtOnce = 1;
184 int rhsCols = b.cols();
187 for(
int k=0; k<rhsCols; k+=NbColsAtOnce)
189 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
190 tmp.
leftCols(actualCols) = b.middleCols(k,actualCols);
198 return *
static_cast<Derived*
>(
this);
202 return *
static_cast<const Derived*
>(
this);
263 template<
typename Rhs>
264 inline const internal::sparse_solve_retval<PastixBase, Rhs>
269 &&
"PastixBase::solve(): invalid number of rows of the right hand side matrix b");
270 return internal::sparse_solve_retval<PastixBase, Rhs>(*
this, b.
derived());
288 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
289 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
314 template <
class Derived>
318 m_iparm.setZero(IPARM_SIZE);
319 m_dparm.setZero(DPARM_SIZE);
321 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
322 pastix(&m_pastixdata, MPI_COMM_WORLD,
324 0, 0, 0, 1, m_iparm.data(), m_dparm.data());
326 m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
327 m_iparm[IPARM_VERBOSE] = 2;
328 m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH;
329 m_iparm[IPARM_INCOMPLETE] = API_NO;
330 m_iparm[IPARM_OOC_LIMIT] = 2000;
331 m_iparm[IPARM_RHS_MAKING] = API_RHS_B;
332 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
334 m_iparm(IPARM_START_TASK) = API_TASK_INIT;
335 m_iparm(IPARM_END_TASK) = API_TASK_INIT;
337 0, 0, 0, 0, m_iparm.data(), m_dparm.data());
340 if(m_iparm(IPARM_ERROR_NUMBER)) {
350 template <
class Derived>
358 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
359 m_isInitialized = m_factorizationIsOk;
363 template <
class Derived>
366 eigen_assert(m_initisOk &&
"The initialization of PaSTiX failed");
373 m_perm.resize(m_size);
374 m_invp.resize(m_size);
376 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
377 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
379 mat.
valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
382 if(m_iparm(IPARM_ERROR_NUMBER))
385 m_analysisIsOk =
false;
390 m_analysisIsOk =
true;
394 template <
class Derived>
398 eigen_assert(m_analysisIsOk &&
"The analysis phase should be called before the factorization phase");
399 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
400 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
404 mat.
valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
407 if(m_iparm(IPARM_ERROR_NUMBER))
410 m_factorizationIsOk =
false;
411 m_isInitialized =
false;
416 m_factorizationIsOk =
true;
417 m_isInitialized =
true;
422 template<
typename Base>
423 template<
typename Rhs,
typename Dest>
426 eigen_assert(m_isInitialized &&
"The matrix should be factorized first");
428 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
433 for (
int i = 0; i < b.cols(); i++){
434 m_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
435 m_iparm[IPARM_END_TASK] = API_TASK_REFINE;
438 m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
444 return m_iparm(IPARM_ERROR_NUMBER)==0;
466 template<
typename _MatrixType,
bool IsStrSym>
473 typedef typename MatrixType::Index
Index;
527 m_iparm(IPARM_SYM) = API_SYM_NO;
528 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
576 template<
typename _MatrixType,
int _UpLo>
631 m_iparm(IPARM_SYM) = API_SYM_YES;
632 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
638 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
657 template<
typename _MatrixType,
int _UpLo>
713 m_iparm(IPARM_SYM) = API_SYM_YES;
714 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
720 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
727 template<
typename _MatrixType,
typename Rhs>
728 struct solve_retval<PastixBase<_MatrixType>, Rhs>
729 : solve_retval_base<PastixBase<_MatrixType>, Rhs>
731 typedef PastixBase<_MatrixType> Dec;
734 template<typename Dest>
void evalTo(Dest& dst)
const
736 dec()._solve(rhs(),dst);
740 template<
typename _MatrixType,
typename Rhs>
741 struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs>
742 : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs>
744 typedef PastixBase<_MatrixType> Dec;
747 template<typename Dest>
void evalTo(Dest& dst)
const
749 dec()._solve_sparse(rhs(),dst);