Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
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 
13 #include "./InternalHeaderCheck.h"
14 
15 namespace Eigen {
16 
17 #if defined(DCOMPLEX)
18  #define PASTIX_COMPLEX COMPLEX
19  #define PASTIX_DCOMPLEX DCOMPLEX
20 #else
21  #define PASTIX_COMPLEX std::complex<float>
22  #define PASTIX_DCOMPLEX std::complex<double>
23 #endif
24 
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;
36 
37 namespace internal
38 {
39 
40  template<class Pastix> struct pastix_traits;
41 
42  template<typename MatrixType_>
43  struct pastix_traits< PastixLU<MatrixType_> >
44  {
45  typedef MatrixType_ MatrixType;
46  typedef typename MatrixType_::Scalar Scalar;
47  typedef typename MatrixType_::RealScalar RealScalar;
48  typedef typename MatrixType_::StorageIndex StorageIndex;
49  };
50 
51  template<typename MatrixType_, int Options>
52  struct pastix_traits< PastixLLT<MatrixType_,Options> >
53  {
54  typedef MatrixType_ MatrixType;
55  typedef typename MatrixType_::Scalar Scalar;
56  typedef typename MatrixType_::RealScalar RealScalar;
57  typedef typename MatrixType_::StorageIndex StorageIndex;
58  };
59 
60  template<typename MatrixType_, int Options>
61  struct pastix_traits< PastixLDLT<MatrixType_,Options> >
62  {
63  typedef MatrixType_ MatrixType;
64  typedef typename MatrixType_::Scalar Scalar;
65  typedef typename MatrixType_::RealScalar RealScalar;
66  typedef typename MatrixType_::StorageIndex StorageIndex;
67  };
68 
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)
70  {
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);
74  }
75 
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)
77  {
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);
81  }
82 
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)
84  {
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);
88  }
89 
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)
91  {
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);
95  }
96 
97  // Convert the matrix to Fortran-style Numbering
98  template <typename MatrixType>
99  void c_to_fortran_numbering (MatrixType& mat)
100  {
101  if ( !(mat.outerIndexPtr()[0]) )
102  {
103  int i;
104  for(i = 0; i <= mat.rows(); ++i)
105  ++mat.outerIndexPtr()[i];
106  for(i = 0; i < mat.nonZeros(); ++i)
107  ++mat.innerIndexPtr()[i];
108  }
109  }
110 
111  // Convert to C-style Numbering
112  template <typename MatrixType>
113  void fortran_to_c_numbering (MatrixType& mat)
114  {
115  // Check the Numbering
116  if ( mat.outerIndexPtr()[0] == 1 )
117  { // Convert to C-style numbering
118  int i;
119  for(i = 0; i <= mat.rows(); ++i)
120  --mat.outerIndexPtr()[i];
121  for(i = 0; i < mat.nonZeros(); ++i)
122  --mat.innerIndexPtr()[i];
123  }
124  }
125 }
126 
127 // This is the base class to interface with PaStiX functions.
128 // Users should not used this class directly.
129 template <class Derived>
130 class PastixBase : public SparseSolverBase<Derived>
131 {
132  protected:
133  typedef SparseSolverBase<Derived> Base;
134  using Base::derived;
135  using Base::m_isInitialized;
136  public:
137  using Base::_solve_impl;
138 
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;
146  enum {
147  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
148  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
149  };
150 
151  public:
152 
153  PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0)
154  {
155  init();
156  }
157 
158  ~PastixBase()
159  {
160  clean();
161  }
162 
163  template<typename Rhs,typename Dest>
164  bool _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
165 
171  Array<StorageIndex,IPARM_SIZE,1>& iparm()
172  {
173  return m_iparm;
174  }
175 
180  int& iparm(int idxparam)
181  {
182  return m_iparm(idxparam);
183  }
184 
189  Array<double,DPARM_SIZE,1>& dparm()
190  {
191  return m_dparm;
192  }
193 
194 
198  double& dparm(int idxparam)
199  {
200  return m_dparm(idxparam);
201  }
202 
203  inline Index cols() const { return m_size; }
204  inline Index rows() const { return m_size; }
205 
214  ComputationInfo info() const
215  {
216  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
217  return m_info;
218  }
219 
220  protected:
221 
222  // Initialize the Pastix data structure, check the matrix
223  void init();
224 
225  // Compute the ordering and the symbolic factorization
226  void analyzePattern(ColSpMatrix& mat);
227 
228  // Compute the numerical factorization
229  void factorize(ColSpMatrix& mat);
230 
231  // Free all the data allocated by Pastix
232  void clean()
233  {
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());
239  }
240 
241  void compute(ColSpMatrix& mat);
242 
243  int m_initisOk;
244  int m_analysisIsOk;
245  int m_factorizationIsOk;
246  mutable ComputationInfo m_info;
247  mutable pastix_data_t *m_pastixdata; // Data structure for pastix
248  mutable int m_comm; // The MPI communicator identifier
249  mutable Array<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
250  mutable Array<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
251  mutable Matrix<StorageIndex,Dynamic,1> m_perm; // Permutation vector
252  mutable Matrix<StorageIndex,Dynamic,1> m_invp; // Inverse permutation vector
253  mutable int m_size; // Size of the matrix
254 };
255 
260 template <class Derived>
261 void PastixBase<Derived>::init()
262 {
263  m_size = 0;
264  m_iparm.setZero(IPARM_SIZE);
265  m_dparm.setZero(DPARM_SIZE);
266 
267  m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
268  pastix(&m_pastixdata, MPI_COMM_WORLD,
269  0, 0, 0, 0,
270  0, 0, 0, 1, m_iparm.data(), m_dparm.data());
271 
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;
279 
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());
284 
285  // Check the returned error
286  if(m_iparm(IPARM_ERROR_NUMBER)) {
287  m_info = InvalidInput;
288  m_initisOk = false;
289  }
290  else {
291  m_info = Success;
292  m_initisOk = true;
293  }
294 }
295 
296 template <class Derived>
297 void PastixBase<Derived>::compute(ColSpMatrix& mat)
298 {
299  eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
300 
301  analyzePattern(mat);
302  factorize(mat);
303 
304  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
305 }
306 
307 
308 template <class Derived>
309 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
310 {
311  eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
312 
313  // clean previous calls
314  if(m_size>0)
315  clean();
316 
317  m_size = internal::convert_index<int>(mat.rows());
318  m_perm.resize(m_size);
319  m_invp.resize(m_size);
320 
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());
325 
326  // Check the returned error
327  if(m_iparm(IPARM_ERROR_NUMBER))
328  {
329  m_info = NumericalIssue;
330  m_analysisIsOk = false;
331  }
332  else
333  {
334  m_info = Success;
335  m_analysisIsOk = true;
336  }
337 }
338 
339 template <class Derived>
340 void PastixBase<Derived>::factorize(ColSpMatrix& mat)
341 {
342 // if(&m_cpyMat != &mat) m_cpyMat = 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());
347 
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());
350 
351  // Check the returned error
352  if(m_iparm(IPARM_ERROR_NUMBER))
353  {
354  m_info = NumericalIssue;
355  m_factorizationIsOk = false;
356  m_isInitialized = false;
357  }
358  else
359  {
360  m_info = Success;
361  m_factorizationIsOk = true;
362  m_isInitialized = true;
363  }
364 }
365 
366 /* Solve the system */
367 template<typename Base>
368 template<typename Rhs,typename Dest>
369 bool PastixBase<Base>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
370 {
371  eigen_assert(m_isInitialized && "The matrix should be factorized first");
372  EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
373  THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
374  int rhs = 1;
375 
376  x = b; /* on return, x is overwritten by the computed solution */
377 
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;
381 
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());
384  }
385 
386  // Check the returned error
387  m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
388 
389  return m_iparm(IPARM_ERROR_NUMBER)==0;
390 }
391 
413 template<typename MatrixType_, bool IsStrSym>
414 class PastixLU : public PastixBase< PastixLU<MatrixType_> >
415 {
416  public:
417  typedef MatrixType_ MatrixType;
418  typedef PastixBase<PastixLU<MatrixType> > Base;
419  typedef typename Base::ColSpMatrix ColSpMatrix;
420  typedef typename MatrixType::StorageIndex StorageIndex;
421 
422  public:
423  PastixLU() : Base()
424  {
425  init();
426  }
427 
428  explicit PastixLU(const MatrixType& matrix):Base()
429  {
430  init();
431  compute(matrix);
432  }
438  void compute (const MatrixType& matrix)
439  {
440  m_structureIsUptodate = false;
441  ColSpMatrix temp;
442  grabMatrix(matrix, temp);
443  Base::compute(temp);
444  }
450  void analyzePattern(const MatrixType& matrix)
451  {
452  m_structureIsUptodate = false;
453  ColSpMatrix temp;
454  grabMatrix(matrix, temp);
455  Base::analyzePattern(temp);
456  }
457 
463  void factorize(const MatrixType& matrix)
464  {
465  ColSpMatrix temp;
466  grabMatrix(matrix, temp);
467  Base::factorize(temp);
468  }
469  protected:
470 
471  void init()
472  {
473  m_structureIsUptodate = false;
474  m_iparm(IPARM_SYM) = API_SYM_NO;
475  m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
476  }
477 
478  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
479  {
480  if(IsStrSym)
481  out = matrix;
482  else
483  {
484  if(!m_structureIsUptodate)
485  {
486  // update the transposed structure
487  m_transposedStructure = matrix.transpose();
488 
489  // Set the elements of the matrix to zero
490  for (Index j=0; j<m_transposedStructure.outerSize(); ++j)
491  for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
492  it.valueRef() = 0.0;
493 
494  m_structureIsUptodate = true;
495  }
496 
497  out = m_transposedStructure + matrix;
498  }
499  internal::c_to_fortran_numbering(out);
500  }
501 
502  using Base::m_iparm;
503  using Base::m_dparm;
504 
505  ColSpMatrix m_transposedStructure;
506  bool m_structureIsUptodate;
507 };
508 
525 template<typename MatrixType_, int UpLo_>
526 class PastixLLT : public PastixBase< PastixLLT<MatrixType_, UpLo_> >
527 {
528  public:
529  typedef MatrixType_ MatrixType;
530  typedef PastixBase<PastixLLT<MatrixType, UpLo_> > Base;
531  typedef typename Base::ColSpMatrix ColSpMatrix;
532 
533  public:
534  enum { UpLo = UpLo_ };
535  PastixLLT() : Base()
536  {
537  init();
538  }
539 
540  explicit PastixLLT(const MatrixType& matrix):Base()
541  {
542  init();
543  compute(matrix);
544  }
545 
549  void compute (const MatrixType& matrix)
550  {
551  ColSpMatrix temp;
552  grabMatrix(matrix, temp);
553  Base::compute(temp);
554  }
555 
560  void analyzePattern(const MatrixType& matrix)
561  {
562  ColSpMatrix temp;
563  grabMatrix(matrix, temp);
564  Base::analyzePattern(temp);
565  }
569  void factorize(const MatrixType& matrix)
570  {
571  ColSpMatrix temp;
572  grabMatrix(matrix, temp);
573  Base::factorize(temp);
574  }
575  protected:
576  using Base::m_iparm;
577 
578  void init()
579  {
580  m_iparm(IPARM_SYM) = API_SYM_YES;
581  m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
582  }
583 
584  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
585  {
586  out.resize(matrix.rows(), matrix.cols());
587  // Pastix supports only lower, column-major matrices
588  out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
589  internal::c_to_fortran_numbering(out);
590  }
591 };
592 
609 template<typename MatrixType_, int UpLo_>
610 class PastixLDLT : public PastixBase< PastixLDLT<MatrixType_, UpLo_> >
611 {
612  public:
613  typedef MatrixType_ MatrixType;
614  typedef PastixBase<PastixLDLT<MatrixType, UpLo_> > Base;
615  typedef typename Base::ColSpMatrix ColSpMatrix;
616 
617  public:
618  enum { UpLo = UpLo_ };
619  PastixLDLT():Base()
620  {
621  init();
622  }
623 
624  explicit PastixLDLT(const MatrixType& matrix):Base()
625  {
626  init();
627  compute(matrix);
628  }
629 
633  void compute (const MatrixType& matrix)
634  {
635  ColSpMatrix temp;
636  grabMatrix(matrix, temp);
637  Base::compute(temp);
638  }
639 
644  void analyzePattern(const MatrixType& matrix)
645  {
646  ColSpMatrix temp;
647  grabMatrix(matrix, temp);
648  Base::analyzePattern(temp);
649  }
653  void factorize(const MatrixType& matrix)
654  {
655  ColSpMatrix temp;
656  grabMatrix(matrix, temp);
657  Base::factorize(temp);
658  }
659 
660  protected:
661  using Base::m_iparm;
662 
663  void init()
664  {
665  m_iparm(IPARM_SYM) = API_SYM_YES;
666  m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
667  }
668 
669  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
670  {
671  // Pastix supports only lower, column-major matrices
672  out.resize(matrix.rows(), matrix.cols());
673  out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
674  internal::c_to_fortran_numbering(out);
675  }
676 };
677 
678 } // end namespace Eigen
679 
680 #endif
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