Eigen  3.2.10
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 #if defined(DCOMPLEX)
14  #define PASTIX_COMPLEX COMPLEX
15  #define PASTIX_DCOMPLEX DCOMPLEX
16 #else
17  #define PASTIX_COMPLEX std::complex<float>
18  #define PASTIX_DCOMPLEX std::complex<double>
19 #endif
20 
21 namespace Eigen {
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::Index Index;
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::Index Index;
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::Index Index;
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 : internal::noncopyable
129 {
130  public:
131  typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
132  typedef _MatrixType MatrixType;
133  typedef typename MatrixType::Scalar Scalar;
134  typedef typename MatrixType::RealScalar RealScalar;
135  typedef typename MatrixType::Index Index;
138 
139  public:
140 
141  PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false), m_pastixdata(0), m_size(0)
142  {
143  init();
144  }
145 
146  ~PastixBase()
147  {
148  clean();
149  }
150 
155  template<typename Rhs>
156  inline const internal::solve_retval<PastixBase, Rhs>
157  solve(const MatrixBase<Rhs>& b) const
158  {
159  eigen_assert(m_isInitialized && "Pastix solver is not initialized.");
160  eigen_assert(rows()==b.rows()
161  && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
162  return internal::solve_retval<PastixBase, Rhs>(*this, b.derived());
163  }
164 
165  template<typename Rhs,typename Dest>
166  bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
167 
168  Derived& derived()
169  {
170  return *static_cast<Derived*>(this);
171  }
172  const Derived& derived() const
173  {
174  return *static_cast<const Derived*>(this);
175  }
176 
183  {
184  return m_iparm;
185  }
186 
191  int& iparm(int idxparam)
192  {
193  return m_iparm(idxparam);
194  }
195 
201  {
202  return m_dparm;
203  }
204 
205 
209  double& dparm(int idxparam)
210  {
211  return m_dparm(idxparam);
212  }
213 
214  inline Index cols() const { return m_size; }
215  inline Index rows() const { return m_size; }
216 
225  ComputationInfo info() const
226  {
227  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
228  return m_info;
229  }
230 
235  template<typename Rhs>
236  inline const internal::sparse_solve_retval<PastixBase, Rhs>
237  solve(const SparseMatrixBase<Rhs>& b) const
238  {
239  eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized.");
240  eigen_assert(rows()==b.rows()
241  && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
242  return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived());
243  }
244 
245  protected:
246 
247  // Initialize the Pastix data structure, check the matrix
248  void init();
249 
250  // Compute the ordering and the symbolic factorization
251  void analyzePattern(ColSpMatrix& mat);
252 
253  // Compute the numerical factorization
254  void factorize(ColSpMatrix& mat);
255 
256  // Free all the data allocated by Pastix
257  void clean()
258  {
259  eigen_assert(m_initisOk && "The Pastix structure should be allocated first");
260  m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
261  m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
262  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
263  m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
264  }
265 
266  void compute(ColSpMatrix& mat);
267 
268  int m_initisOk;
269  int m_analysisIsOk;
270  int m_factorizationIsOk;
271  bool m_isInitialized;
272  mutable ComputationInfo m_info;
273  mutable pastix_data_t *m_pastixdata; // Data structure for pastix
274  mutable int m_comm; // The MPI communicator identifier
275  mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
276  mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
277  mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector
278  mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector
279  mutable int m_size; // Size of the matrix
280 };
281 
286 template <class Derived>
287 void PastixBase<Derived>::init()
288 {
289  m_size = 0;
290  m_iparm.setZero(IPARM_SIZE);
291  m_dparm.setZero(DPARM_SIZE);
292 
293  m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
294  pastix(&m_pastixdata, MPI_COMM_WORLD,
295  0, 0, 0, 0,
296  0, 0, 0, 1, m_iparm.data(), m_dparm.data());
297 
298  m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
299  m_iparm[IPARM_VERBOSE] = 2;
300  m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH;
301  m_iparm[IPARM_INCOMPLETE] = API_NO;
302  m_iparm[IPARM_OOC_LIMIT] = 2000;
303  m_iparm[IPARM_RHS_MAKING] = API_RHS_B;
304  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
305 
306  m_iparm(IPARM_START_TASK) = API_TASK_INIT;
307  m_iparm(IPARM_END_TASK) = API_TASK_INIT;
308  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
309  0, 0, 0, 0, m_iparm.data(), m_dparm.data());
310 
311  // Check the returned error
312  if(m_iparm(IPARM_ERROR_NUMBER)) {
313  m_info = InvalidInput;
314  m_initisOk = false;
315  }
316  else {
317  m_info = Success;
318  m_initisOk = true;
319  }
320 }
321 
322 template <class Derived>
323 void PastixBase<Derived>::compute(ColSpMatrix& mat)
324 {
325  eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
326 
327  analyzePattern(mat);
328  factorize(mat);
329 
330  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
331  m_isInitialized = m_factorizationIsOk;
332 }
333 
334 
335 template <class Derived>
336 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
337 {
338  eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
339 
340  // clean previous calls
341  if(m_size>0)
342  clean();
343 
344  m_size = mat.rows();
345  m_perm.resize(m_size);
346  m_invp.resize(m_size);
347 
348  m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
349  m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
350  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
351  mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
352 
353  // Check the returned error
354  if(m_iparm(IPARM_ERROR_NUMBER))
355  {
356  m_info = NumericalIssue;
357  m_analysisIsOk = false;
358  }
359  else
360  {
361  m_info = Success;
362  m_analysisIsOk = true;
363  }
364 }
365 
366 template <class Derived>
367 void PastixBase<Derived>::factorize(ColSpMatrix& mat)
368 {
369 // if(&m_cpyMat != &mat) m_cpyMat = mat;
370  eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
371  m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
372  m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
373  m_size = mat.rows();
374 
375  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
376  mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
377 
378  // Check the returned error
379  if(m_iparm(IPARM_ERROR_NUMBER))
380  {
381  m_info = NumericalIssue;
382  m_factorizationIsOk = false;
383  m_isInitialized = false;
384  }
385  else
386  {
387  m_info = Success;
388  m_factorizationIsOk = true;
389  m_isInitialized = true;
390  }
391 }
392 
393 /* Solve the system */
394 template<typename Base>
395 template<typename Rhs,typename Dest>
396 bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
397 {
398  eigen_assert(m_isInitialized && "The matrix should be factorized first");
399  EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
400  THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
401  int rhs = 1;
402 
403  x = b; /* on return, x is overwritten by the computed solution */
404 
405  for (int i = 0; i < b.cols(); i++){
406  m_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
407  m_iparm[IPARM_END_TASK] = API_TASK_REFINE;
408 
409  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0,
410  m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
411  }
412 
413  // Check the returned error
414  m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
415 
416  return m_iparm(IPARM_ERROR_NUMBER)==0;
417 }
418 
438 template<typename _MatrixType, bool IsStrSym>
439 class PastixLU : public PastixBase< PastixLU<_MatrixType> >
440 {
441  public:
442  typedef _MatrixType MatrixType;
443  typedef PastixBase<PastixLU<MatrixType> > Base;
444  typedef typename Base::ColSpMatrix ColSpMatrix;
445  typedef typename MatrixType::Index Index;
446 
447  public:
448  PastixLU() : Base()
449  {
450  init();
451  }
452 
453  PastixLU(const MatrixType& matrix):Base()
454  {
455  init();
456  compute(matrix);
457  }
463  void compute (const MatrixType& matrix)
464  {
465  m_structureIsUptodate = false;
466  ColSpMatrix temp;
467  grabMatrix(matrix, temp);
468  Base::compute(temp);
469  }
475  void analyzePattern(const MatrixType& matrix)
476  {
477  m_structureIsUptodate = false;
478  ColSpMatrix temp;
479  grabMatrix(matrix, temp);
480  Base::analyzePattern(temp);
481  }
482 
488  void factorize(const MatrixType& matrix)
489  {
490  ColSpMatrix temp;
491  grabMatrix(matrix, temp);
492  Base::factorize(temp);
493  }
494  protected:
495 
496  void init()
497  {
498  m_structureIsUptodate = false;
499  m_iparm(IPARM_SYM) = API_SYM_NO;
500  m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
501  }
502 
503  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
504  {
505  if(IsStrSym)
506  out = matrix;
507  else
508  {
509  if(!m_structureIsUptodate)
510  {
511  // update the transposed structure
512  m_transposedStructure = matrix.transpose();
513 
514  // Set the elements of the matrix to zero
515  for (Index j=0; j<m_transposedStructure.outerSize(); ++j)
516  for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
517  it.valueRef() = 0.0;
518 
519  m_structureIsUptodate = true;
520  }
521 
522  out = m_transposedStructure + matrix;
523  }
524  internal::c_to_fortran_numbering(out);
525  }
526 
527  using Base::m_iparm;
528  using Base::m_dparm;
529 
530  ColSpMatrix m_transposedStructure;
531  bool m_structureIsUptodate;
532 };
533 
548 template<typename _MatrixType, int _UpLo>
549 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
550 {
551  public:
552  typedef _MatrixType MatrixType;
553  typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
554  typedef typename Base::ColSpMatrix ColSpMatrix;
555 
556  public:
557  enum { UpLo = _UpLo };
558  PastixLLT() : Base()
559  {
560  init();
561  }
562 
563  PastixLLT(const MatrixType& matrix):Base()
564  {
565  init();
566  compute(matrix);
567  }
568 
572  void compute (const MatrixType& matrix)
573  {
574  ColSpMatrix temp;
575  grabMatrix(matrix, temp);
576  Base::compute(temp);
577  }
578 
583  void analyzePattern(const MatrixType& matrix)
584  {
585  ColSpMatrix temp;
586  grabMatrix(matrix, temp);
587  Base::analyzePattern(temp);
588  }
592  void factorize(const MatrixType& matrix)
593  {
594  ColSpMatrix temp;
595  grabMatrix(matrix, temp);
596  Base::factorize(temp);
597  }
598  protected:
599  using Base::m_iparm;
600 
601  void init()
602  {
603  m_iparm(IPARM_SYM) = API_SYM_YES;
604  m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
605  }
606 
607  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
608  {
609  // Pastix supports only lower, column-major matrices
610  out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
611  internal::c_to_fortran_numbering(out);
612  }
613 };
614 
629 template<typename _MatrixType, int _UpLo>
630 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
631 {
632  public:
633  typedef _MatrixType MatrixType;
634  typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base;
635  typedef typename Base::ColSpMatrix ColSpMatrix;
636 
637  public:
638  enum { UpLo = _UpLo };
639  PastixLDLT():Base()
640  {
641  init();
642  }
643 
644  PastixLDLT(const MatrixType& matrix):Base()
645  {
646  init();
647  compute(matrix);
648  }
649 
653  void compute (const MatrixType& matrix)
654  {
655  ColSpMatrix temp;
656  grabMatrix(matrix, temp);
657  Base::compute(temp);
658  }
659 
664  void analyzePattern(const MatrixType& matrix)
665  {
666  ColSpMatrix temp;
667  grabMatrix(matrix, temp);
668  Base::analyzePattern(temp);
669  }
673  void factorize(const MatrixType& matrix)
674  {
675  ColSpMatrix temp;
676  grabMatrix(matrix, temp);
677  Base::factorize(temp);
678  }
679 
680  protected:
681  using Base::m_iparm;
682 
683  void init()
684  {
685  m_iparm(IPARM_SYM) = API_SYM_YES;
686  m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
687  }
688 
689  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
690  {
691  // Pastix supports only lower, column-major matrices
692  out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
693  internal::c_to_fortran_numbering(out);
694  }
695 };
696 
697 namespace internal {
698 
699 template<typename _MatrixType, typename Rhs>
700 struct solve_retval<PastixBase<_MatrixType>, Rhs>
701  : solve_retval_base<PastixBase<_MatrixType>, Rhs>
702 {
703  typedef PastixBase<_MatrixType> Dec;
704  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
705 
706  template<typename Dest> void evalTo(Dest& dst) const
707  {
708  dec()._solve(rhs(),dst);
709  }
710 };
711 
712 template<typename _MatrixType, typename Rhs>
713 struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs>
714  : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs>
715 {
716  typedef PastixBase<_MatrixType> Dec;
717  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
718 
719  template<typename Dest> void evalTo(Dest& dst) const
720  {
721  this->defaultEvalTo(dst);
722  }
723 };
724 
725 } // end namespace internal
726 
727 } // end namespace Eigen
728 
729 #endif
A versatible sparse matrix representation.
Definition: SparseMatrix.h:85
Definition: Constants.h:378
Definition: LDLT.h:16
const Scalar * valuePtr() const
Definition: SparseMatrix.h:131
Index rows() const
Definition: SparseMatrixBase.h:160
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:239
A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library...
Definition: PaStiXSupport.h:33
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:664
Derived & derived()
Definition: EigenBase.h:34
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:583
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:653
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:673
Definition: Constants.h:383
Index rows() const
Definition: SparseMatrix.h:119
Definition: Eigen_Colamd.h:50
void analyzePattern(const MatrixType &matrix)
Definition: PaStiXSupport.h:475
const Index * outerIndexPtr() const
Definition: SparseMatrix.h:149
Index cols() const
Definition: SparseMatrix.h:121
const Index * innerIndexPtr() const
Definition: SparseMatrix.h:140
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:42
Definition: Constants.h:376
const unsigned int RowMajorBit
Definition: Constants.h:53
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
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:592
ComputationInfo
Definition: Constants.h:374
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:463
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
void factorize(const MatrixType &matrix)
Definition: PaStiXSupport.h:488
void compute(const MatrixType &matrix)
Definition: PaStiXSupport.h:572
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515