Eigen  3.3.7
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  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.
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  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 
258 template <class Derived>
259 void 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 
294 template <class Derived>
295 void 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 
306 template <class Derived>
307 void 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 
337 template <class Derived>
338 void 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 */
365 template<typename Base>
366 template<typename Rhs,typename Dest>
367 bool 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 
411 template<typename _MatrixType, bool IsStrSym>
412 class 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 
523 template<typename _MatrixType, int _UpLo>
524 class 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 
607 template<typename _MatrixType, int _UpLo>
608 class 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
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:436
A versatible sparse matrix representation.
Definition: SparseMatrix.h:96
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:461
Namespace containing all symbols from the Eigen library.
Definition: Core:306
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:651
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:547
const unsigned int RowMajorBit
Definition: Constants.h:61
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:631
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:558
A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library...
Definition: PaStiXSupport.h:33
Definition: Constants.h:434
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Definition: Constants.h:439
Definition: Constants.h:432
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:448
Definition: Eigen_Colamd.h:50
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:567
Index outerSize() const
Definition: SparseMatrix.h:143
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:430
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:642