10 #ifndef EIGEN_PASTIXSUPPORT_H
11 #define EIGEN_PASTIXSUPPORT_H
16 #define PASTIX_COMPLEX COMPLEX
17 #define PASTIX_DCOMPLEX DCOMPLEX
19 #define PASTIX_COMPLEX std::complex<float>
20 #define PASTIX_DCOMPLEX std::complex<double>
31 template<
typename _MatrixType,
bool IsStrSym = false>
class PastixLU;
32 template<
typename _MatrixType,
int Options>
class PastixLLT;
33 template<
typename _MatrixType,
int Options>
class PastixLDLT;
38 template<
class Pastix>
struct pastix_traits;
40 template<
typename _MatrixType>
41 struct pastix_traits<
PastixLU<_MatrixType> >
43 typedef _MatrixType MatrixType;
44 typedef typename _MatrixType::Scalar Scalar;
45 typedef typename _MatrixType::RealScalar RealScalar;
46 typedef typename _MatrixType::StorageIndex StorageIndex;
49 template<
typename _MatrixType,
int Options>
50 struct pastix_traits<
PastixLLT<_MatrixType,Options> >
52 typedef _MatrixType MatrixType;
53 typedef typename _MatrixType::Scalar Scalar;
54 typedef typename _MatrixType::RealScalar RealScalar;
55 typedef typename _MatrixType::StorageIndex StorageIndex;
58 template<
typename _MatrixType,
int Options>
59 struct pastix_traits< PastixLDLT<_MatrixType,Options> >
61 typedef _MatrixType MatrixType;
62 typedef typename _MatrixType::Scalar Scalar;
63 typedef typename _MatrixType::RealScalar RealScalar;
64 typedef typename _MatrixType::StorageIndex StorageIndex;
67 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)
69 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
70 if (nbrhs == 0) {x = NULL; nbrhs=1;}
71 s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
74 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)
76 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
77 if (nbrhs == 0) {x = NULL; nbrhs=1;}
78 d_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, std::complex<float> *vals,
int *perm,
int * invp, std::complex<float> *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 c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm);
88 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)
90 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
91 if (nbrhs == 0) {x = NULL; nbrhs=1;}
92 z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm);
96 template <
typename MatrixType>
97 void c_to_fortran_numbering (MatrixType& mat)
99 if ( !(mat.outerIndexPtr()[0]) )
102 for(i = 0; i <= mat.rows(); ++i)
103 ++mat.outerIndexPtr()[i];
104 for(i = 0; i < mat.nonZeros(); ++i)
105 ++mat.innerIndexPtr()[i];
110 template <
typename MatrixType>
111 void fortran_to_c_numbering (MatrixType& mat)
114 if ( mat.outerIndexPtr()[0] == 1 )
117 for(i = 0; i <= mat.rows(); ++i)
118 --mat.outerIndexPtr()[i];
119 for(i = 0; i < mat.nonZeros(); ++i)
120 --mat.innerIndexPtr()[i];
127 template <
class Derived>
128 class PastixBase :
public SparseSolverBase<Derived>
131 typedef SparseSolverBase<Derived> Base;
133 using Base::m_isInitialized;
135 using Base::_solve_impl;
137 typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
138 typedef _MatrixType MatrixType;
139 typedef typename MatrixType::Scalar Scalar;
140 typedef typename MatrixType::RealScalar RealScalar;
141 typedef typename MatrixType::StorageIndex StorageIndex;
142 typedef Matrix<Scalar,Dynamic,1> Vector;
143 typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix;
147 PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0)
157 template<
typename Rhs,
typename Dest>
158 bool _solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &x)
const;
165 Array<StorageIndex,IPARM_SIZE,1>& iparm()
174 int& iparm(
int idxparam)
176 return m_iparm(idxparam);
183 Array<RealScalar,IPARM_SIZE,1>& dparm()
192 double& dparm(
int idxparam)
194 return m_dparm(idxparam);
197 inline Index cols()
const {
return m_size; }
198 inline Index rows()
const {
return m_size; }
210 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
220 void analyzePattern(ColSpMatrix& mat);
223 void factorize(ColSpMatrix& mat);
228 eigen_assert(m_initisOk &&
"The Pastix structure should be allocated first");
229 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
230 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
231 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
232 m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
235 void compute(ColSpMatrix& mat);
239 int m_factorizationIsOk;
241 mutable pastix_data_t *m_pastixdata;
243 mutable Matrix<int,IPARM_SIZE,1> m_iparm;
244 mutable Matrix<double,DPARM_SIZE,1> m_dparm;
245 mutable Matrix<StorageIndex,Dynamic,1> m_perm;
246 mutable Matrix<StorageIndex,Dynamic,1> m_invp;
254 template <
class Derived>
255 void PastixBase<Derived>::init()
258 m_iparm.setZero(IPARM_SIZE);
259 m_dparm.setZero(DPARM_SIZE);
261 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
262 pastix(&m_pastixdata, MPI_COMM_WORLD,
264 0, 0, 0, 1, m_iparm.data(), m_dparm.data());
266 m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
267 m_iparm[IPARM_VERBOSE] = 2;
268 m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH;
269 m_iparm[IPARM_INCOMPLETE] = API_NO;
270 m_iparm[IPARM_OOC_LIMIT] = 2000;
271 m_iparm[IPARM_RHS_MAKING] = API_RHS_B;
272 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
274 m_iparm(IPARM_START_TASK) = API_TASK_INIT;
275 m_iparm(IPARM_END_TASK) = API_TASK_INIT;
276 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
277 0, 0, 0, 0, m_iparm.data(), m_dparm.data());
280 if(m_iparm(IPARM_ERROR_NUMBER)) {
290 template <
class Derived>
291 void PastixBase<Derived>::compute(ColSpMatrix& mat)
293 eigen_assert(mat.rows() == mat.cols() &&
"The input matrix should be squared");
298 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
302 template <
class Derived>
303 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
305 eigen_assert(m_initisOk &&
"The initialization of PaSTiX failed");
311 m_size = internal::convert_index<int>(mat.rows());
312 m_perm.resize(m_size);
313 m_invp.resize(m_size);
315 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
316 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
317 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
318 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
321 if(m_iparm(IPARM_ERROR_NUMBER))
324 m_analysisIsOk =
false;
329 m_analysisIsOk =
true;
333 template <
class Derived>
334 void PastixBase<Derived>::factorize(ColSpMatrix& mat)
337 eigen_assert(m_analysisIsOk &&
"The analysis phase should be called before the factorization phase");
338 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
339 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
340 m_size = internal::convert_index<int>(mat.rows());
342 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
343 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
346 if(m_iparm(IPARM_ERROR_NUMBER))
349 m_factorizationIsOk =
false;
350 m_isInitialized =
false;
355 m_factorizationIsOk =
true;
356 m_isInitialized =
true;
361 template<
typename Base>
362 template<
typename Rhs,
typename Dest>
363 bool PastixBase<Base>::_solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &x)
const
365 eigen_assert(m_isInitialized &&
"The matrix should be factorized first");
367 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
372 for (
int i = 0; i < b.cols(); i++){
373 m_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
374 m_iparm[IPARM_END_TASK] = API_TASK_REFINE;
376 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, internal::convert_index<int>(x.rows()), 0, 0, 0,
377 m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
383 return m_iparm(IPARM_ERROR_NUMBER)==0;
405 template<
typename _MatrixType,
bool IsStrSym>
406 class PastixLU :
public PastixBase< PastixLU<_MatrixType> >
409 typedef _MatrixType MatrixType;
410 typedef PastixBase<PastixLU<MatrixType> > Base;
411 typedef typename Base::ColSpMatrix ColSpMatrix;
412 typedef typename MatrixType::StorageIndex StorageIndex;
420 explicit PastixLU(
const MatrixType& matrix):Base()
432 m_structureIsUptodate =
false;
434 grabMatrix(matrix, temp);
444 m_structureIsUptodate =
false;
446 grabMatrix(matrix, temp);
447 Base::analyzePattern(temp);
458 grabMatrix(matrix, temp);
459 Base::factorize(temp);
465 m_structureIsUptodate =
false;
466 m_iparm(IPARM_SYM) = API_SYM_NO;
467 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
470 void grabMatrix(
const MatrixType& matrix, ColSpMatrix& out)
476 if(!m_structureIsUptodate)
479 m_transposedStructure = matrix.transpose();
482 for (Index j=0; j<m_transposedStructure.
outerSize(); ++j)
483 for(
typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
486 m_structureIsUptodate =
true;
489 out = m_transposedStructure + matrix;
491 internal::c_to_fortran_numbering(out);
497 ColSpMatrix m_transposedStructure;
498 bool m_structureIsUptodate;
515 template<
typename _MatrixType,
int _UpLo>
516 class PastixLLT :
public PastixBase< PastixLLT<_MatrixType, _UpLo> >
519 typedef _MatrixType MatrixType;
520 typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
521 typedef typename Base::ColSpMatrix ColSpMatrix;
524 enum { UpLo = _UpLo };
530 explicit PastixLLT(
const MatrixType& matrix):Base()
542 grabMatrix(matrix, temp);
553 grabMatrix(matrix, temp);
554 Base::analyzePattern(temp);
562 grabMatrix(matrix, temp);
563 Base::factorize(temp);
570 m_iparm(IPARM_SYM) = API_SYM_YES;
571 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
574 void grabMatrix(
const MatrixType& matrix, ColSpMatrix& out)
577 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
578 internal::c_to_fortran_numbering(out);
596 template<
typename _MatrixType,
int _UpLo>
597 class PastixLDLT :
public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
600 typedef _MatrixType MatrixType;
601 typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base;
602 typedef typename Base::ColSpMatrix ColSpMatrix;
605 enum { UpLo = _UpLo };
611 explicit PastixLDLT(
const MatrixType& matrix):Base()
623 grabMatrix(matrix, temp);
634 grabMatrix(matrix, temp);
635 Base::analyzePattern(temp);
643 grabMatrix(matrix, temp);
644 Base::factorize(temp);
652 m_iparm(IPARM_SYM) = API_SYM_YES;
653 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
656 void grabMatrix(
const MatrixType& matrix, ColSpMatrix& out)
659 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
660 internal::c_to_fortran_numbering(out);
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:430
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:455
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:640
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:539
Index outerSize() const
Definition: SparseMatrix.h:137
const unsigned int RowMajorBit
Definition: Constants.h:53
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:620
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:550
A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library...
Definition: PaStiXSupport.h:33
Definition: Constants.h:426
Definition: Constants.h:431
Definition: Constants.h:424
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:442
Definition: Eigen_Colamd.h:54
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:559
Interface to the PaStix solver.
Definition: PaStiXSupport.h:31
A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library...
Definition: PaStiXSupport.h:32
ComputationInfo
Definition: Constants.h:422
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:631