Eigen  3.2.91
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 namespace 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 
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;
34 
35 namespace 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  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  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  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  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.
127 template <class Derived>
128 class 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 
145  public:
146 
147  PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0)
148  {
149  init();
150  }
151 
152  ~PastixBase()
153  {
154  clean();
155  }
156 
157  template<typename Rhs,typename Dest>
158  bool _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
159 
165  Array<StorageIndex,IPARM_SIZE,1>& iparm()
166  {
167  return m_iparm;
168  }
169 
174  int& iparm(int idxparam)
175  {
176  return m_iparm(idxparam);
177  }
178 
183  Array<RealScalar,IPARM_SIZE,1>& dparm()
184  {
185  return m_dparm;
186  }
187 
188 
192  double& dparm(int idxparam)
193  {
194  return m_dparm(idxparam);
195  }
196 
197  inline Index cols() const { return m_size; }
198  inline Index rows() const { return m_size; }
199 
208  ComputationInfo info() const
209  {
210  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
211  return m_info;
212  }
213 
214  protected:
215 
216  // Initialize the Pastix data structure, check the matrix
217  void init();
218 
219  // Compute the ordering and the symbolic factorization
220  void analyzePattern(ColSpMatrix& mat);
221 
222  // Compute the numerical factorization
223  void factorize(ColSpMatrix& mat);
224 
225  // Free all the data allocated by Pastix
226  void clean()
227  {
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());
233  }
234 
235  void compute(ColSpMatrix& mat);
236 
237  int m_initisOk;
238  int m_analysisIsOk;
239  int m_factorizationIsOk;
240  mutable ComputationInfo m_info;
241  mutable pastix_data_t *m_pastixdata; // Data structure for pastix
242  mutable int m_comm; // The MPI communicator identifier
243  mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
244  mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
245  mutable Matrix<StorageIndex,Dynamic,1> m_perm; // Permutation vector
246  mutable Matrix<StorageIndex,Dynamic,1> m_invp; // Inverse permutation vector
247  mutable int m_size; // Size of the matrix
248 };
249 
254 template <class Derived>
255 void PastixBase<Derived>::init()
256 {
257  m_size = 0;
258  m_iparm.setZero(IPARM_SIZE);
259  m_dparm.setZero(DPARM_SIZE);
260 
261  m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
262  pastix(&m_pastixdata, MPI_COMM_WORLD,
263  0, 0, 0, 0,
264  0, 0, 0, 1, m_iparm.data(), m_dparm.data());
265 
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;
273 
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());
278 
279  // Check the returned error
280  if(m_iparm(IPARM_ERROR_NUMBER)) {
281  m_info = InvalidInput;
282  m_initisOk = false;
283  }
284  else {
285  m_info = Success;
286  m_initisOk = true;
287  }
288 }
289 
290 template <class Derived>
291 void PastixBase<Derived>::compute(ColSpMatrix& mat)
292 {
293  eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
294 
295  analyzePattern(mat);
296  factorize(mat);
297 
298  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
299 }
300 
301 
302 template <class Derived>
303 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
304 {
305  eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
306 
307  // clean previous calls
308  if(m_size>0)
309  clean();
310 
311  m_size = internal::convert_index<int>(mat.rows());
312  m_perm.resize(m_size);
313  m_invp.resize(m_size);
314 
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());
319 
320  // Check the returned error
321  if(m_iparm(IPARM_ERROR_NUMBER))
322  {
323  m_info = NumericalIssue;
324  m_analysisIsOk = false;
325  }
326  else
327  {
328  m_info = Success;
329  m_analysisIsOk = true;
330  }
331 }
332 
333 template <class Derived>
334 void PastixBase<Derived>::factorize(ColSpMatrix& mat)
335 {
336 // if(&m_cpyMat != &mat) m_cpyMat = 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());
341 
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());
344 
345  // Check the returned error
346  if(m_iparm(IPARM_ERROR_NUMBER))
347  {
348  m_info = NumericalIssue;
349  m_factorizationIsOk = false;
350  m_isInitialized = false;
351  }
352  else
353  {
354  m_info = Success;
355  m_factorizationIsOk = true;
356  m_isInitialized = true;
357  }
358 }
359 
360 /* Solve the system */
361 template<typename Base>
362 template<typename Rhs,typename Dest>
363 bool PastixBase<Base>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
364 {
365  eigen_assert(m_isInitialized && "The matrix should be factorized first");
366  EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
367  THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
368  int rhs = 1;
369 
370  x = b; /* on return, x is overwritten by the computed solution */
371 
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;
375 
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());
378  }
379 
380  // Check the returned error
381  m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
382 
383  return m_iparm(IPARM_ERROR_NUMBER)==0;
384 }
385 
405 template<typename _MatrixType, bool IsStrSym>
406 class PastixLU : public PastixBase< PastixLU<_MatrixType> >
407 {
408  public:
409  typedef _MatrixType MatrixType;
410  typedef PastixBase<PastixLU<MatrixType> > Base;
411  typedef typename Base::ColSpMatrix ColSpMatrix;
412  typedef typename MatrixType::StorageIndex StorageIndex;
413 
414  public:
415  PastixLU() : Base()
416  {
417  init();
418  }
419 
420  explicit PastixLU(const MatrixType& matrix):Base()
421  {
422  init();
423  compute(matrix);
424  }
430  void compute (const MatrixType& matrix)
431  {
432  m_structureIsUptodate = false;
433  ColSpMatrix temp;
434  grabMatrix(matrix, temp);
435  Base::compute(temp);
436  }
442  void analyzePattern(const MatrixType& matrix)
443  {
444  m_structureIsUptodate = false;
445  ColSpMatrix temp;
446  grabMatrix(matrix, temp);
447  Base::analyzePattern(temp);
448  }
449 
455  void factorize(const MatrixType& matrix)
456  {
457  ColSpMatrix temp;
458  grabMatrix(matrix, temp);
459  Base::factorize(temp);
460  }
461  protected:
462 
463  void init()
464  {
465  m_structureIsUptodate = false;
466  m_iparm(IPARM_SYM) = API_SYM_NO;
467  m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
468  }
469 
470  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
471  {
472  if(IsStrSym)
473  out = matrix;
474  else
475  {
476  if(!m_structureIsUptodate)
477  {
478  // update the transposed structure
479  m_transposedStructure = matrix.transpose();
480 
481  // Set the elements of the matrix to zero
482  for (Index j=0; j<m_transposedStructure.outerSize(); ++j)
483  for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
484  it.valueRef() = 0.0;
485 
486  m_structureIsUptodate = true;
487  }
488 
489  out = m_transposedStructure + matrix;
490  }
491  internal::c_to_fortran_numbering(out);
492  }
493 
494  using Base::m_iparm;
495  using Base::m_dparm;
496 
497  ColSpMatrix m_transposedStructure;
498  bool m_structureIsUptodate;
499 };
500 
515 template<typename _MatrixType, int _UpLo>
516 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
517 {
518  public:
519  typedef _MatrixType MatrixType;
520  typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
521  typedef typename Base::ColSpMatrix ColSpMatrix;
522 
523  public:
524  enum { UpLo = _UpLo };
525  PastixLLT() : Base()
526  {
527  init();
528  }
529 
530  explicit PastixLLT(const MatrixType& matrix):Base()
531  {
532  init();
533  compute(matrix);
534  }
535 
539  void compute (const MatrixType& matrix)
540  {
541  ColSpMatrix temp;
542  grabMatrix(matrix, temp);
543  Base::compute(temp);
544  }
545 
550  void analyzePattern(const MatrixType& matrix)
551  {
552  ColSpMatrix temp;
553  grabMatrix(matrix, temp);
554  Base::analyzePattern(temp);
555  }
559  void factorize(const MatrixType& matrix)
560  {
561  ColSpMatrix temp;
562  grabMatrix(matrix, temp);
563  Base::factorize(temp);
564  }
565  protected:
566  using Base::m_iparm;
567 
568  void init()
569  {
570  m_iparm(IPARM_SYM) = API_SYM_YES;
571  m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
572  }
573 
574  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
575  {
576  // Pastix supports only lower, column-major matrices
577  out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
578  internal::c_to_fortran_numbering(out);
579  }
580 };
581 
596 template<typename _MatrixType, int _UpLo>
597 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
598 {
599  public:
600  typedef _MatrixType MatrixType;
601  typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base;
602  typedef typename Base::ColSpMatrix ColSpMatrix;
603 
604  public:
605  enum { UpLo = _UpLo };
606  PastixLDLT():Base()
607  {
608  init();
609  }
610 
611  explicit PastixLDLT(const MatrixType& matrix):Base()
612  {
613  init();
614  compute(matrix);
615  }
616 
620  void compute (const MatrixType& matrix)
621  {
622  ColSpMatrix temp;
623  grabMatrix(matrix, temp);
624  Base::compute(temp);
625  }
626 
631  void analyzePattern(const MatrixType& matrix)
632  {
633  ColSpMatrix temp;
634  grabMatrix(matrix, temp);
635  Base::analyzePattern(temp);
636  }
640  void factorize(const MatrixType& matrix)
641  {
642  ColSpMatrix temp;
643  grabMatrix(matrix, temp);
644  Base::factorize(temp);
645  }
646 
647  protected:
648  using Base::m_iparm;
649 
650  void init()
651  {
652  m_iparm(IPARM_SYM) = API_SYM_YES;
653  m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
654  }
655 
656  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
657  {
658  // Pastix supports only lower, column-major matrices
659  out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
660  internal::c_to_fortran_numbering(out);
661  }
662 };
663 
664 } // end namespace Eigen
665 
666 #endif
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:430
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:455
Definition: LDLT.h:16
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