PaStiXSupport.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_PASTIXSUPPORT_H
26 #define EIGEN_PASTIXSUPPORT_H
27 
28 namespace Eigen {
29 
38 template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
39 template<typename _MatrixType, int Options> class PastixLLT;
40 template<typename _MatrixType, int Options> class PastixLDLT;
41 
42 namespace internal
43 {
44 
45  template<class Pastix> struct pastix_traits;
46 
47  template<typename _MatrixType>
48  struct pastix_traits< PastixLU<_MatrixType> >
49  {
50  typedef _MatrixType MatrixType;
51  typedef typename _MatrixType::Scalar Scalar;
52  typedef typename _MatrixType::RealScalar RealScalar;
53  typedef typename _MatrixType::Index Index;
54  };
55 
56  template<typename _MatrixType, int Options>
57  struct pastix_traits< PastixLLT<_MatrixType,Options> >
58  {
59  typedef _MatrixType MatrixType;
60  typedef typename _MatrixType::Scalar Scalar;
61  typedef typename _MatrixType::RealScalar RealScalar;
62  typedef typename _MatrixType::Index Index;
63  };
64 
65  template<typename _MatrixType, int Options>
66  struct pastix_traits< PastixLDLT<_MatrixType,Options> >
67  {
68  typedef _MatrixType MatrixType;
69  typedef typename _MatrixType::Scalar Scalar;
70  typedef typename _MatrixType::RealScalar RealScalar;
71  typedef typename _MatrixType::Index Index;
72  };
73 
74  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)
75  {
76  if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
77  if (nbrhs == 0) {x = NULL; nbrhs=1;}
78  s_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, double *vals, int *perm, int * invp, double *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  d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, 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<float> *vals, int *perm, int * invp, std::complex<float> *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  c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm);
93  }
94 
95  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)
96  {
97  if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
98  if (nbrhs == 0) {x = NULL; nbrhs=1;}
99  z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm);
100  }
101 
102  // Convert the matrix to Fortran-style Numbering
103  template <typename MatrixType>
104  void c_to_fortran_numbering (MatrixType& mat)
105  {
106  if ( !(mat.outerIndexPtr()[0]) )
107  {
108  int i;
109  for(i = 0; i <= mat.rows(); ++i)
110  ++mat.outerIndexPtr()[i];
111  for(i = 0; i < mat.nonZeros(); ++i)
112  ++mat.innerIndexPtr()[i];
113  }
114  }
115 
116  // Convert to C-style Numbering
117  template <typename MatrixType>
118  void fortran_to_c_numbering (MatrixType& mat)
119  {
120  // Check the Numbering
121  if ( mat.outerIndexPtr()[0] == 1 )
122  { // Convert to C-style numbering
123  int i;
124  for(i = 0; i <= mat.rows(); ++i)
125  --mat.outerIndexPtr()[i];
126  for(i = 0; i < mat.nonZeros(); ++i)
127  --mat.innerIndexPtr()[i];
128  }
129  }
130 }
131 
132 // This is the base class to interface with PaStiX functions.
133 // Users should not used this class directly.
134 template <class Derived>
135 class PastixBase : internal::noncopyable
136 {
137  public:
138  typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
140  typedef typename MatrixType::Scalar Scalar;
141  typedef typename MatrixType::RealScalar RealScalar;
142  typedef typename MatrixType::Index Index;
145 
146  public:
147 
149  {
150  init();
151  }
152 
154  {
155  clean();
156  }
157 
162  template<typename Rhs>
163  inline const internal::solve_retval<PastixBase, Rhs>
164  solve(const MatrixBase<Rhs>& b) const
165  {
166  eigen_assert(m_isInitialized && "Pastix solver is not initialized.");
167  eigen_assert(rows()==b.rows()
168  && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
169  return internal::solve_retval<PastixBase, Rhs>(*this, b.derived());
170  }
171 
172  template<typename Rhs,typename Dest>
173  bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
174 
176  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
178  {
179  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
180  eigen_assert(rows()==b.rows());
181 
182  // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
183  static const int NbColsAtOnce = 1;
184  int rhsCols = b.cols();
185  int size = b.rows();
187  for(int k=0; k<rhsCols; k+=NbColsAtOnce)
188  {
189  int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
190  tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
191  tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
192  dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
193  }
194  }
195 
196  Derived& derived()
197  {
198  return *static_cast<Derived*>(this);
199  }
200  const Derived& derived() const
201  {
202  return *static_cast<const Derived*>(this);
203  }
204 
211  {
212  return m_iparm;
213  }
214 
219  int& iparm(int idxparam)
220  {
221  return m_iparm(idxparam);
222  }
223 
229  {
230  return m_dparm;
231  }
232 
233 
237  double& dparm(int idxparam)
238  {
239  return m_dparm(idxparam);
240  }
241 
242  inline Index cols() const { return m_size; }
243  inline Index rows() const { return m_size; }
244 
254  {
255  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
256  return m_info;
257  }
258 
263  template<typename Rhs>
264  inline const internal::sparse_solve_retval<PastixBase, Rhs>
265  solve(const SparseMatrixBase<Rhs>& b) const
266  {
267  eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized.");
268  eigen_assert(rows()==b.rows()
269  && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
270  return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived());
271  }
272 
273  protected:
274 
275  // Initialize the Pastix data structure, check the matrix
276  void init();
277 
278  // Compute the ordering and the symbolic factorization
279  void analyzePattern(ColSpMatrix& mat);
280 
281  // Compute the numerical factorization
282  void factorize(ColSpMatrix& mat);
283 
284  // Free all the data allocated by Pastix
285  void clean()
286  {
287  eigen_assert(m_initisOk && "The Pastix structure should be allocated first");
288  m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
289  m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
290  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
291  m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
292  }
293 
294  void compute(ColSpMatrix& mat);
295 
301  mutable pastix_data_t *m_pastixdata; // Data structure for pastix
302  mutable int m_comm; // The MPI communicator identifier
303  mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
304  mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
305  mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector
306  mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector
307  mutable int m_size; // Size of the matrix
308 };
309 
314 template <class Derived>
316 {
317  m_size = 0;
318  m_iparm.setZero(IPARM_SIZE);
319  m_dparm.setZero(DPARM_SIZE);
320 
321  m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
322  pastix(&m_pastixdata, MPI_COMM_WORLD,
323  0, 0, 0, 0,
324  0, 0, 0, 1, m_iparm.data(), m_dparm.data());
325 
326  m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
327  m_iparm[IPARM_VERBOSE] = 2;
328  m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH;
329  m_iparm[IPARM_INCOMPLETE] = API_NO;
330  m_iparm[IPARM_OOC_LIMIT] = 2000;
331  m_iparm[IPARM_RHS_MAKING] = API_RHS_B;
332  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
333 
334  m_iparm(IPARM_START_TASK) = API_TASK_INIT;
335  m_iparm(IPARM_END_TASK) = API_TASK_INIT;
336  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
337  0, 0, 0, 0, m_iparm.data(), m_dparm.data());
338 
339  // Check the returned error
340  if(m_iparm(IPARM_ERROR_NUMBER)) {
341  m_info = InvalidInput;
342  m_initisOk = false;
343  }
344  else {
345  m_info = Success;
346  m_initisOk = true;
347  }
348 }
349 
350 template <class Derived>
352 {
353  eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
354 
355  analyzePattern(mat);
356  factorize(mat);
357 
358  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
359  m_isInitialized = m_factorizationIsOk;
360 }
361 
362 
363 template <class Derived>
365 {
366  eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
367 
368  // clean previous calls
369  if(m_size>0)
370  clean();
371 
372  m_size = mat.rows();
373  m_perm.resize(m_size);
374  m_invp.resize(m_size);
375 
376  m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
377  m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
378  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
379  mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
380 
381  // Check the returned error
382  if(m_iparm(IPARM_ERROR_NUMBER))
383  {
384  m_info = NumericalIssue;
385  m_analysisIsOk = false;
386  }
387  else
388  {
389  m_info = Success;
390  m_analysisIsOk = true;
391  }
392 }
393 
394 template <class Derived>
396 {
397 // if(&m_cpyMat != &mat) m_cpyMat = mat;
398  eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
399  m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
400  m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
401  m_size = mat.rows();
402 
403  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
404  mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
405 
406  // Check the returned error
407  if(m_iparm(IPARM_ERROR_NUMBER))
408  {
409  m_info = NumericalIssue;
410  m_factorizationIsOk = false;
411  m_isInitialized = false;
412  }
413  else
414  {
415  m_info = Success;
416  m_factorizationIsOk = true;
417  m_isInitialized = true;
418  }
419 }
420 
421 /* Solve the system */
422 template<typename Base>
423 template<typename Rhs,typename Dest>
425 {
426  eigen_assert(m_isInitialized && "The matrix should be factorized first");
427  EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
428  THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
429  int rhs = 1;
430 
431  x = b; /* on return, x is overwritten by the computed solution */
432 
433  for (int i = 0; i < b.cols(); i++){
434  m_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
435  m_iparm[IPARM_END_TASK] = API_TASK_REFINE;
436 
437  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0,
438  m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
439  }
440 
441  // Check the returned error
442  m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
443 
444  return m_iparm(IPARM_ERROR_NUMBER)==0;
445 }
446 
466 template<typename _MatrixType, bool IsStrSym>
467 class PastixLU : public PastixBase< PastixLU<_MatrixType> >
468 {
469  public:
472  typedef typename Base::ColSpMatrix ColSpMatrix;
473  typedef typename MatrixType::Index Index;
474 
475  public:
477  {
478  init();
479  }
480 
481  PastixLU(const MatrixType& matrix):Base()
482  {
483  init();
484  compute(matrix);
485  }
491  void compute (const MatrixType& matrix)
492  {
493  m_structureIsUptodate = false;
494  ColSpMatrix temp;
495  grabMatrix(matrix, temp);
496  Base::compute(temp);
497  }
503  void analyzePattern(const MatrixType& matrix)
504  {
505  m_structureIsUptodate = false;
506  ColSpMatrix temp;
507  grabMatrix(matrix, temp);
508  Base::analyzePattern(temp);
509  }
510 
516  void factorize(const MatrixType& matrix)
517  {
518  ColSpMatrix temp;
519  grabMatrix(matrix, temp);
520  Base::factorize(temp);
521  }
522  protected:
523 
524  void init()
525  {
526  m_structureIsUptodate = false;
527  m_iparm(IPARM_SYM) = API_SYM_NO;
528  m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
529  }
530 
531  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
532  {
533  if(IsStrSym)
534  out = matrix;
535  else
536  {
538  {
539  // update the transposed structure
540  m_transposedStructure = matrix.transpose();
541 
542  // Set the elements of the matrix to zero
543  for (Index j=0; j<m_transposedStructure.outerSize(); ++j)
544  for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
545  it.valueRef() = 0.0;
546 
547  m_structureIsUptodate = true;
548  }
549 
550  out = m_transposedStructure + matrix;
551  }
553  }
554 
555  using Base::m_iparm;
556  using Base::m_dparm;
557 
560 };
561 
576 template<typename _MatrixType, int _UpLo>
577 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
578 {
579  public:
582  typedef typename Base::ColSpMatrix ColSpMatrix;
583 
584  public:
585  enum { UpLo = _UpLo };
587  {
588  init();
589  }
590 
591  PastixLLT(const MatrixType& matrix):Base()
592  {
593  init();
594  compute(matrix);
595  }
596 
600  void compute (const MatrixType& matrix)
601  {
602  ColSpMatrix temp;
603  grabMatrix(matrix, temp);
604  Base::compute(temp);
605  }
606 
611  void analyzePattern(const MatrixType& matrix)
612  {
613  ColSpMatrix temp;
614  grabMatrix(matrix, temp);
615  Base::analyzePattern(temp);
616  }
620  void factorize(const MatrixType& matrix)
621  {
622  ColSpMatrix temp;
623  grabMatrix(matrix, temp);
624  Base::factorize(temp);
625  }
626  protected:
627  using Base::m_iparm;
628 
629  void init()
630  {
631  m_iparm(IPARM_SYM) = API_SYM_YES;
632  m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
633  }
634 
635  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
636  {
637  // Pastix supports only lower, column-major matrices
638  out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
640  }
641 };
642 
657 template<typename _MatrixType, int _UpLo>
658 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
659 {
660  public:
663  typedef typename Base::ColSpMatrix ColSpMatrix;
664 
665  public:
666  enum { UpLo = _UpLo };
668  {
669  init();
670  }
671 
672  PastixLDLT(const MatrixType& matrix):Base()
673  {
674  init();
675  compute(matrix);
676  }
677 
681  void compute (const MatrixType& matrix)
682  {
683  ColSpMatrix temp;
684  grabMatrix(matrix, temp);
685  Base::compute(temp);
686  }
687 
692  void analyzePattern(const MatrixType& matrix)
693  {
694  ColSpMatrix temp;
695  grabMatrix(matrix, temp);
696  Base::analyzePattern(temp);
697  }
701  void factorize(const MatrixType& matrix)
702  {
703  ColSpMatrix temp;
704  grabMatrix(matrix, temp);
705  Base::factorize(temp);
706  }
707 
708  protected:
709  using Base::m_iparm;
710 
711  void init()
712  {
713  m_iparm(IPARM_SYM) = API_SYM_YES;
714  m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
715  }
716 
717  void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
718  {
719  // Pastix supports only lower, column-major matrices
720  out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
722  }
723 };
724 
725 namespace internal {
726 
727 template<typename _MatrixType, typename Rhs>
728 struct solve_retval<PastixBase<_MatrixType>, Rhs>
729  : solve_retval_base<PastixBase<_MatrixType>, Rhs>
730 {
731  typedef PastixBase<_MatrixType> Dec;
732  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
733 
734  template<typename Dest> void evalTo(Dest& dst) const
735  {
736  dec()._solve(rhs(),dst);
737  }
738 };
739 
740 template<typename _MatrixType, typename Rhs>
741 struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs>
742  : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs>
743 {
744  typedef PastixBase<_MatrixType> Dec;
746 
747  template<typename Dest> void evalTo(Dest& dst) const
748  {
749  dec()._solve_sparse(rhs(),dst);
750  }
751 };
752 
753 } // end namespace internal
754 
755 } // end namespace Eigen
756 
757 #endif