Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
PaStiXSupport.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_PASTIXSUPPORT_H
11#define EIGEN_PASTIXSUPPORT_H
12
13namespace Eigen {
14
15#if defined(DCOMPLEX)
16 #define PASTIX_COMPLEX COMPLEX
17 #define PASTIX_DCOMPLEX DCOMPLEX
18#else
19 #define PASTIX_COMPLEX std::complex<float>
20 #define PASTIX_DCOMPLEX std::complex<double>
21#endif
22
31template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
32template<typename _MatrixType, int Options> class PastixLLT;
33template<typename _MatrixType, int Options> class PastixLDLT;
34
35namespace internal
36{
37
38 template<class Pastix> struct pastix_traits;
39
40 template<typename _MatrixType>
41 struct pastix_traits< PastixLU<_MatrixType> >
42 {
43 typedef _MatrixType MatrixType;
44 typedef typename _MatrixType::Scalar Scalar;
45 typedef typename _MatrixType::RealScalar RealScalar;
46 typedef typename _MatrixType::StorageIndex StorageIndex;
47 };
48
49 template<typename _MatrixType, int Options>
50 struct pastix_traits< PastixLLT<_MatrixType,Options> >
51 {
52 typedef _MatrixType MatrixType;
53 typedef typename _MatrixType::Scalar Scalar;
54 typedef typename _MatrixType::RealScalar RealScalar;
55 typedef typename _MatrixType::StorageIndex StorageIndex;
56 };
57
58 template<typename _MatrixType, int Options>
59 struct pastix_traits< PastixLDLT<_MatrixType,Options> >
60 {
61 typedef _MatrixType MatrixType;
62 typedef typename _MatrixType::Scalar Scalar;
63 typedef typename _MatrixType::RealScalar RealScalar;
64 typedef typename _MatrixType::StorageIndex StorageIndex;
65 };
66
67 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)
68 {
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);
72 }
73
74 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)
75 {
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);
79 }
80
81 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)
82 {
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);
86 }
87
88 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)
89 {
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);
93 }
94
95 // Convert the matrix to Fortran-style Numbering
96 template <typename MatrixType>
97 void c_to_fortran_numbering (MatrixType& mat)
98 {
99 if ( !(mat.outerIndexPtr()[0]) )
100 {
101 int i;
102 for(i = 0; i <= mat.rows(); ++i)
103 ++mat.outerIndexPtr()[i];
104 for(i = 0; i < mat.nonZeros(); ++i)
105 ++mat.innerIndexPtr()[i];
106 }
107 }
108
109 // Convert to C-style Numbering
110 template <typename MatrixType>
111 void fortran_to_c_numbering (MatrixType& mat)
112 {
113 // Check the Numbering
114 if ( mat.outerIndexPtr()[0] == 1 )
115 { // Convert to C-style numbering
116 int i;
117 for(i = 0; i <= mat.rows(); ++i)
118 --mat.outerIndexPtr()[i];
119 for(i = 0; i < mat.nonZeros(); ++i)
120 --mat.innerIndexPtr()[i];
121 }
122 }
123}
124
125// This is the base class to interface with PaStiX functions.
126// Users should not used this class directly.
127template <class Derived>
128class PastixBase : public SparseSolverBase<Derived>
129{
130 protected:
131 typedef SparseSolverBase<Derived> Base;
132 using Base::derived;
133 using Base::m_isInitialized;
134 public:
135 using Base::_solve_impl;
136
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;
144 enum {
145 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
146 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
147 };
148
149 public:
150
151 PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0)
152 {
153 init();
154 }
155
156 ~PastixBase()
157 {
158 clean();
159 }
160
161 template<typename Rhs,typename Dest>
162 bool _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
163
169 Array<StorageIndex,IPARM_SIZE,1>& iparm()
170 {
171 return m_iparm;
172 }
173
178 int& iparm(int idxparam)
179 {
180 return m_iparm(idxparam);
181 }
182
187 Array<double,DPARM_SIZE,1>& dparm()
188 {
189 return m_dparm;
190 }
191
192
196 double& dparm(int idxparam)
197 {
198 return m_dparm(idxparam);
199 }
200
201 inline Index cols() const { return m_size; }
202 inline Index rows() const { return m_size; }
203
212 ComputationInfo info() const
213 {
214 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
215 return m_info;
216 }
217
218 protected:
219
220 // Initialize the Pastix data structure, check the matrix
221 void init();
222
223 // Compute the ordering and the symbolic factorization
224 void analyzePattern(ColSpMatrix& mat);
225
226 // Compute the numerical factorization
227 void factorize(ColSpMatrix& mat);
228
229 // Free all the data allocated by Pastix
230 void clean()
231 {
232 eigen_assert(m_initisOk && "The Pastix structure should be allocated first");
233 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
234 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
235 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
236 m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
237 }
238
239 void compute(ColSpMatrix& mat);
240
241 int m_initisOk;
242 int m_analysisIsOk;
243 int m_factorizationIsOk;
244 mutable ComputationInfo m_info;
245 mutable pastix_data_t *m_pastixdata; // Data structure for pastix
246 mutable int m_comm; // The MPI communicator identifier
247 mutable Array<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
248 mutable Array<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
249 mutable Matrix<StorageIndex,Dynamic,1> m_perm; // Permutation vector
250 mutable Matrix<StorageIndex,Dynamic,1> m_invp; // Inverse permutation vector
251 mutable int m_size; // Size of the matrix
252};
253
258template <class Derived>
259void PastixBase<Derived>::init()
260{
261 m_size = 0;
262 m_iparm.setZero(IPARM_SIZE);
263 m_dparm.setZero(DPARM_SIZE);
264
265 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
266 pastix(&m_pastixdata, MPI_COMM_WORLD,
267 0, 0, 0, 0,
268 0, 0, 0, 1, m_iparm.data(), m_dparm.data());
269
270 m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
271 m_iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
272 m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH;
273 m_iparm[IPARM_INCOMPLETE] = API_NO;
274 m_iparm[IPARM_OOC_LIMIT] = 2000;
275 m_iparm[IPARM_RHS_MAKING] = API_RHS_B;
276 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
277
278 m_iparm(IPARM_START_TASK) = API_TASK_INIT;
279 m_iparm(IPARM_END_TASK) = API_TASK_INIT;
280 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
281 0, 0, 0, 0, m_iparm.data(), m_dparm.data());
282
283 // Check the returned error
284 if(m_iparm(IPARM_ERROR_NUMBER)) {
285 m_info = InvalidInput;
286 m_initisOk = false;
287 }
288 else {
289 m_info = Success;
290 m_initisOk = true;
291 }
292}
293
294template <class Derived>
295void PastixBase<Derived>::compute(ColSpMatrix& mat)
296{
297 eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
298
299 analyzePattern(mat);
300 factorize(mat);
301
302 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
303}
304
305
306template <class Derived>
307void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
308{
309 eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
310
311 // clean previous calls
312 if(m_size>0)
313 clean();
314
315 m_size = internal::convert_index<int>(mat.rows());
316 m_perm.resize(m_size);
317 m_invp.resize(m_size);
318
319 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
320 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
321 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
322 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
323
324 // Check the returned error
325 if(m_iparm(IPARM_ERROR_NUMBER))
326 {
327 m_info = NumericalIssue;
328 m_analysisIsOk = false;
329 }
330 else
331 {
332 m_info = Success;
333 m_analysisIsOk = true;
334 }
335}
336
337template <class Derived>
338void PastixBase<Derived>::factorize(ColSpMatrix& mat)
339{
340// if(&m_cpyMat != &mat) m_cpyMat = mat;
341 eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
342 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
343 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
344 m_size = internal::convert_index<int>(mat.rows());
345
346 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
347 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
348
349 // Check the returned error
350 if(m_iparm(IPARM_ERROR_NUMBER))
351 {
352 m_info = NumericalIssue;
353 m_factorizationIsOk = false;
354 m_isInitialized = false;
355 }
356 else
357 {
358 m_info = Success;
359 m_factorizationIsOk = true;
360 m_isInitialized = true;
361 }
362}
363
364/* Solve the system */
365template<typename Base>
366template<typename Rhs,typename Dest>
367bool PastixBase<Base>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
368{
369 eigen_assert(m_isInitialized && "The matrix should be factorized first");
370 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
371 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
372 int rhs = 1;
373
374 x = b; /* on return, x is overwritten by the computed solution */
375
376 for (int i = 0; i < b.cols(); i++){
377 m_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
378 m_iparm[IPARM_END_TASK] = API_TASK_REFINE;
379
380 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, internal::convert_index<int>(x.rows()), 0, 0, 0,
381 m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
382 }
383
384 // Check the returned error
385 m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
386
387 return m_iparm(IPARM_ERROR_NUMBER)==0;
388}
389
411template<typename _MatrixType, bool IsStrSym>
412class PastixLU : public PastixBase< PastixLU<_MatrixType> >
413{
414 public:
415 typedef _MatrixType MatrixType;
416 typedef PastixBase<PastixLU<MatrixType> > Base;
417 typedef typename Base::ColSpMatrix ColSpMatrix;
418 typedef typename MatrixType::StorageIndex StorageIndex;
419
420 public:
421 PastixLU() : Base()
422 {
423 init();
424 }
425
426 explicit PastixLU(const MatrixType& matrix):Base()
427 {
428 init();
429 compute(matrix);
430 }
436 void compute (const MatrixType& matrix)
437 {
438 m_structureIsUptodate = false;
439 ColSpMatrix temp;
440 grabMatrix(matrix, temp);
441 Base::compute(temp);
442 }
448 void analyzePattern(const MatrixType& matrix)
449 {
450 m_structureIsUptodate = false;
451 ColSpMatrix temp;
452 grabMatrix(matrix, temp);
453 Base::analyzePattern(temp);
454 }
455
461 void factorize(const MatrixType& matrix)
462 {
463 ColSpMatrix temp;
464 grabMatrix(matrix, temp);
465 Base::factorize(temp);
466 }
467 protected:
468
469 void init()
470 {
471 m_structureIsUptodate = false;
472 m_iparm(IPARM_SYM) = API_SYM_NO;
473 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
474 }
475
476 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
477 {
478 if(IsStrSym)
479 out = matrix;
480 else
481 {
482 if(!m_structureIsUptodate)
483 {
484 // update the transposed structure
485 m_transposedStructure = matrix.transpose();
486
487 // Set the elements of the matrix to zero
488 for (Index j=0; j<m_transposedStructure.outerSize(); ++j)
489 for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
490 it.valueRef() = 0.0;
491
492 m_structureIsUptodate = true;
493 }
494
495 out = m_transposedStructure + matrix;
496 }
497 internal::c_to_fortran_numbering(out);
498 }
499
500 using Base::m_iparm;
501 using Base::m_dparm;
502
503 ColSpMatrix m_transposedStructure;
504 bool m_structureIsUptodate;
505};
506
523template<typename _MatrixType, int _UpLo>
524class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
525{
526 public:
527 typedef _MatrixType MatrixType;
528 typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
529 typedef typename Base::ColSpMatrix ColSpMatrix;
530
531 public:
532 enum { UpLo = _UpLo };
533 PastixLLT() : Base()
534 {
535 init();
536 }
537
538 explicit PastixLLT(const MatrixType& matrix):Base()
539 {
540 init();
541 compute(matrix);
542 }
543
547 void compute (const MatrixType& matrix)
548 {
549 ColSpMatrix temp;
550 grabMatrix(matrix, temp);
551 Base::compute(temp);
552 }
553
558 void analyzePattern(const MatrixType& matrix)
559 {
560 ColSpMatrix temp;
561 grabMatrix(matrix, temp);
562 Base::analyzePattern(temp);
563 }
567 void factorize(const MatrixType& matrix)
568 {
569 ColSpMatrix temp;
570 grabMatrix(matrix, temp);
571 Base::factorize(temp);
572 }
573 protected:
574 using Base::m_iparm;
575
576 void init()
577 {
578 m_iparm(IPARM_SYM) = API_SYM_YES;
579 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
580 }
581
582 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
583 {
584 out.resize(matrix.rows(), matrix.cols());
585 // Pastix supports only lower, column-major matrices
586 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
587 internal::c_to_fortran_numbering(out);
588 }
589};
590
607template<typename _MatrixType, int _UpLo>
608class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
609{
610 public:
611 typedef _MatrixType MatrixType;
612 typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base;
613 typedef typename Base::ColSpMatrix ColSpMatrix;
614
615 public:
616 enum { UpLo = _UpLo };
617 PastixLDLT():Base()
618 {
619 init();
620 }
621
622 explicit PastixLDLT(const MatrixType& matrix):Base()
623 {
624 init();
625 compute(matrix);
626 }
627
631 void compute (const MatrixType& matrix)
632 {
633 ColSpMatrix temp;
634 grabMatrix(matrix, temp);
635 Base::compute(temp);
636 }
637
642 void analyzePattern(const MatrixType& matrix)
643 {
644 ColSpMatrix temp;
645 grabMatrix(matrix, temp);
646 Base::analyzePattern(temp);
647 }
651 void factorize(const MatrixType& matrix)
652 {
653 ColSpMatrix temp;
654 grabMatrix(matrix, temp);
655 Base::factorize(temp);
656 }
657
658 protected:
659 using Base::m_iparm;
660
661 void init()
662 {
663 m_iparm(IPARM_SYM) = API_SYM_YES;
664 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
665 }
666
667 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
668 {
669 // Pastix supports only lower, column-major matrices
670 out.resize(matrix.rows(), matrix.cols());
671 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
672 internal::c_to_fortran_numbering(out);
673 }
674};
675
676} // end namespace Eigen
677
678#endif
A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library.
Definition: PaStiXSupport.h:609
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:642
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:651
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:631
A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library.
Definition: PaStiXSupport.h:525
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:547
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:567
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:558
Interface to the PaStix solver.
Definition: PaStiXSupport.h:413
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:448
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:461
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:436
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
Index outerSize() const
Definition: SparseMatrix.h:145
ComputationInfo
Definition: Constants.h:440
@ NumericalIssue
Definition: Constants.h:444
@ InvalidInput
Definition: Constants.h:449
@ Success
Definition: Constants.h:442
Matrix< Type, Size, 1 > Vector
[c++11]
Definition: Matrix.h:551
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74