10 #ifndef EIGEN_PASTIXSUPPORT_H
11 #define EIGEN_PASTIXSUPPORT_H
13 #include "./InternalHeaderCheck.h"
18 #define PASTIX_COMPLEX COMPLEX
19 #define PASTIX_DCOMPLEX DCOMPLEX
21 #define PASTIX_COMPLEX std::complex<float>
22 #define PASTIX_DCOMPLEX std::complex<double>
33 template<
typename MatrixType_,
bool IsStrSym = false>
class PastixLU;
34 template<
typename MatrixType_,
int Options>
class PastixLLT;
35 template<
typename MatrixType_,
int Options>
class PastixLDLT;
40 template<
class Pastix>
struct pastix_traits;
42 template<
typename MatrixType_>
43 struct pastix_traits< PastixLU<MatrixType_> >
45 typedef MatrixType_ MatrixType;
46 typedef typename MatrixType_::Scalar Scalar;
47 typedef typename MatrixType_::RealScalar RealScalar;
48 typedef typename MatrixType_::StorageIndex StorageIndex;
51 template<
typename MatrixType_,
int Options>
52 struct pastix_traits< PastixLLT<MatrixType_,Options> >
54 typedef MatrixType_ MatrixType;
55 typedef typename MatrixType_::Scalar Scalar;
56 typedef typename MatrixType_::RealScalar RealScalar;
57 typedef typename MatrixType_::StorageIndex StorageIndex;
60 template<
typename MatrixType_,
int Options>
61 struct pastix_traits< PastixLDLT<MatrixType_,Options> >
63 typedef MatrixType_ MatrixType;
64 typedef typename MatrixType_::Scalar Scalar;
65 typedef typename MatrixType_::RealScalar RealScalar;
66 typedef typename MatrixType_::StorageIndex StorageIndex;
69 inline 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)
71 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
72 if (nbrhs == 0) {x = NULL; nbrhs=1;}
73 s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
76 inline 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)
78 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
79 if (nbrhs == 0) {x = NULL; nbrhs=1;}
80 d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
83 inline 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)
85 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
86 if (nbrhs == 0) {x = NULL; nbrhs=1;}
87 c_pastix(pastix_data, pastix_comm, n, ptr, idx,
reinterpret_cast<PASTIX_COMPLEX*
>(vals), perm, invp,
reinterpret_cast<PASTIX_COMPLEX*
>(x), nbrhs, iparm, dparm);
90 inline 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)
92 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
93 if (nbrhs == 0) {x = NULL; nbrhs=1;}
94 z_pastix(pastix_data, pastix_comm, n, ptr, idx,
reinterpret_cast<PASTIX_DCOMPLEX*
>(vals), perm, invp,
reinterpret_cast<PASTIX_DCOMPLEX*
>(x), nbrhs, iparm, dparm);
98 template <
typename MatrixType>
99 void c_to_fortran_numbering (MatrixType& mat)
101 if ( !(mat.outerIndexPtr()[0]) )
104 for(i = 0; i <= mat.rows(); ++i)
105 ++mat.outerIndexPtr()[i];
106 for(i = 0; i < mat.nonZeros(); ++i)
107 ++mat.innerIndexPtr()[i];
112 template <
typename MatrixType>
113 void fortran_to_c_numbering (MatrixType& mat)
116 if ( mat.outerIndexPtr()[0] == 1 )
119 for(i = 0; i <= mat.rows(); ++i)
120 --mat.outerIndexPtr()[i];
121 for(i = 0; i < mat.nonZeros(); ++i)
122 --mat.innerIndexPtr()[i];
129 template <
class Derived>
130 class PastixBase :
public SparseSolverBase<Derived>
133 typedef SparseSolverBase<Derived> Base;
135 using Base::m_isInitialized;
137 using Base::_solve_impl;
139 typedef typename internal::pastix_traits<Derived>::MatrixType MatrixType_;
140 typedef MatrixType_ MatrixType;
141 typedef typename MatrixType::Scalar Scalar;
142 typedef typename MatrixType::RealScalar RealScalar;
143 typedef typename MatrixType::StorageIndex StorageIndex;
144 typedef Matrix<Scalar,Dynamic,1>
Vector;
145 typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix;
147 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
148 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
153 PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0)
163 template<
typename Rhs,
typename Dest>
164 bool _solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &x)
const;
171 Array<StorageIndex,IPARM_SIZE,1>& iparm()
180 int& iparm(
int idxparam)
182 return m_iparm(idxparam);
189 Array<double,DPARM_SIZE,1>& dparm()
198 double& dparm(
int idxparam)
200 return m_dparm(idxparam);
203 inline Index cols()
const {
return m_size; }
204 inline Index rows()
const {
return m_size; }
216 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
226 void analyzePattern(ColSpMatrix& mat);
229 void factorize(ColSpMatrix& mat);
234 eigen_assert(m_initisOk &&
"The Pastix structure should be allocated first");
235 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
236 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
237 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
238 m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
241 void compute(ColSpMatrix& mat);
245 int m_factorizationIsOk;
247 mutable pastix_data_t *m_pastixdata;
249 mutable Array<int,IPARM_SIZE,1> m_iparm;
250 mutable Array<double,DPARM_SIZE,1> m_dparm;
251 mutable Matrix<StorageIndex,Dynamic,1> m_perm;
252 mutable Matrix<StorageIndex,Dynamic,1> m_invp;
260 template <
class Derived>
261 void PastixBase<Derived>::init()
264 m_iparm.setZero(IPARM_SIZE);
265 m_dparm.setZero(DPARM_SIZE);
267 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
268 pastix(&m_pastixdata, MPI_COMM_WORLD,
270 0, 0, 0, 1, m_iparm.data(), m_dparm.data());
272 m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
273 m_iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
274 m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH;
275 m_iparm[IPARM_INCOMPLETE] = API_NO;
276 m_iparm[IPARM_OOC_LIMIT] = 2000;
277 m_iparm[IPARM_RHS_MAKING] = API_RHS_B;
278 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
280 m_iparm(IPARM_START_TASK) = API_TASK_INIT;
281 m_iparm(IPARM_END_TASK) = API_TASK_INIT;
282 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
283 0, 0, 0, 0, m_iparm.data(), m_dparm.data());
286 if(m_iparm(IPARM_ERROR_NUMBER)) {
296 template <
class Derived>
297 void PastixBase<Derived>::compute(ColSpMatrix& mat)
299 eigen_assert(mat.rows() == mat.cols() &&
"The input matrix should be squared");
304 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
308 template <
class Derived>
309 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
311 eigen_assert(m_initisOk &&
"The initialization of PaSTiX failed");
317 m_size = internal::convert_index<int>(mat.rows());
318 m_perm.resize(m_size);
319 m_invp.resize(m_size);
321 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
322 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
323 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
324 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
327 if(m_iparm(IPARM_ERROR_NUMBER))
330 m_analysisIsOk =
false;
335 m_analysisIsOk =
true;
339 template <
class Derived>
340 void PastixBase<Derived>::factorize(ColSpMatrix& mat)
343 eigen_assert(m_analysisIsOk &&
"The analysis phase should be called before the factorization phase");
344 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
345 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
346 m_size = internal::convert_index<int>(mat.rows());
348 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
349 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
352 if(m_iparm(IPARM_ERROR_NUMBER))
355 m_factorizationIsOk =
false;
356 m_isInitialized =
false;
361 m_factorizationIsOk =
true;
362 m_isInitialized =
true;
367 template<
typename Base>
368 template<
typename Rhs,
typename Dest>
369 bool PastixBase<Base>::_solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &x)
const
371 eigen_assert(m_isInitialized &&
"The matrix should be factorized first");
373 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
378 for (
int i = 0; i < b.cols(); i++){
379 m_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
380 m_iparm[IPARM_END_TASK] = API_TASK_REFINE;
382 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, internal::convert_index<int>(x.rows()), 0, 0, 0,
383 m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
389 return m_iparm(IPARM_ERROR_NUMBER)==0;
413 template<
typename MatrixType_,
bool IsStrSym>
414 class PastixLU :
public PastixBase< PastixLU<MatrixType_> >
417 typedef MatrixType_ MatrixType;
418 typedef PastixBase<PastixLU<MatrixType> >
Base;
420 typedef typename MatrixType::StorageIndex StorageIndex;
440 m_structureIsUptodate =
false;
442 grabMatrix(matrix, temp);
452 m_structureIsUptodate =
false;
454 grabMatrix(matrix, temp);
455 Base::analyzePattern(temp);
466 grabMatrix(matrix, temp);
467 Base::factorize(temp);
473 m_structureIsUptodate =
false;
474 m_iparm(IPARM_SYM) = API_SYM_NO;
475 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
478 void grabMatrix(
const MatrixType& matrix, ColSpMatrix& out)
484 if(!m_structureIsUptodate)
487 m_transposedStructure = matrix.transpose();
491 for(
typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
494 m_structureIsUptodate =
true;
497 out = m_transposedStructure + matrix;
499 internal::c_to_fortran_numbering(out);
505 ColSpMatrix m_transposedStructure;
506 bool m_structureIsUptodate;
525 template<
typename MatrixType_,
int UpLo_>
526 class PastixLLT :
public PastixBase< PastixLLT<MatrixType_, UpLo_> >
529 typedef MatrixType_ MatrixType;
530 typedef PastixBase<PastixLLT<MatrixType, UpLo_> >
Base;
534 enum { UpLo = UpLo_ };
552 grabMatrix(matrix, temp);
563 grabMatrix(matrix, temp);
564 Base::analyzePattern(temp);
572 grabMatrix(matrix, temp);
573 Base::factorize(temp);
580 m_iparm(IPARM_SYM) = API_SYM_YES;
581 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
584 void grabMatrix(
const MatrixType& matrix, ColSpMatrix& out)
586 out.resize(matrix.rows(), matrix.cols());
588 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
589 internal::c_to_fortran_numbering(out);
609 template<
typename MatrixType_,
int UpLo_>
610 class PastixLDLT :
public PastixBase< PastixLDLT<MatrixType_, UpLo_> >
613 typedef MatrixType_ MatrixType;
614 typedef PastixBase<PastixLDLT<MatrixType, UpLo_> >
Base;
618 enum { UpLo = UpLo_ };
636 grabMatrix(matrix, temp);
647 grabMatrix(matrix, temp);
648 Base::analyzePattern(temp);
656 grabMatrix(matrix, temp);
657 Base::factorize(temp);
665 m_iparm(IPARM_SYM) = API_SYM_YES;
666 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
669 void grabMatrix(
const MatrixType& matrix, ColSpMatrix& out)
672 out.resize(matrix.rows(), matrix.cols());
673 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
674 internal::c_to_fortran_numbering(out);
A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library.
Definition: PaStiXSupport.h:611
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:633
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:644
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:653
A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library.
Definition: PaStiXSupport.h:527
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:549
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:569
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:560
Interface to the PaStix solver.
Definition: PaStiXSupport.h:415
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:450
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:463
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:438
A versatible sparse matrix representation.
Definition: SparseMatrix.h:100
Index outerSize() const
Definition: SparseMatrix.h:147
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ NumericalIssue
Definition: Constants.h:446
@ InvalidInput
Definition: Constants.h:451
@ Success
Definition: Constants.h:444
const unsigned int RowMajorBit
Definition: Constants.h:68
Matrix< Type, Size, 1 > Vector
[c++11] Size×1 vector of type Type.
Definition: Matrix.h:544
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59