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 
39 template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
40 template<typename _MatrixType, int Options> class PastixLLT;
41 template<typename _MatrixType, int Options> class PastixLDLT;
42 
43 namespace internal
44 {
45 
46  template<class Pastix> struct pastix_traits;
47 
48  template<typename _MatrixType>
49  struct pastix_traits< PastixLU<_MatrixType> >
50  {
51  typedef _MatrixType MatrixType;
52  typedef typename _MatrixType::Scalar Scalar;
53  typedef typename _MatrixType::RealScalar RealScalar;
54  typedef typename _MatrixType::Index Index;
55  };
56 
57  template<typename _MatrixType, int Options>
58  struct pastix_traits< PastixLLT<_MatrixType,Options> >
59  {
60  typedef _MatrixType MatrixType;
61  typedef typename _MatrixType::Scalar Scalar;
62  typedef typename _MatrixType::RealScalar RealScalar;
63  typedef typename _MatrixType::Index Index;
64  };
65 
66  template<typename _MatrixType, int Options>
67  struct pastix_traits< PastixLDLT<_MatrixType,Options> >
68  {
69  typedef _MatrixType MatrixType;
70  typedef typename _MatrixType::Scalar Scalar;
71  typedef typename _MatrixType::RealScalar RealScalar;
72  typedef typename _MatrixType::Index Index;
73  };
74 
75  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)
76  {
77  if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
78  if (nbrhs == 0) x = NULL;
79  s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
80  }
81 
82  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)
83  {
84  if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
85  if (nbrhs == 0) x = NULL;
86  d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
87  }
88 
89  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)
90  {
91  c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm);
92  }
93 
94  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)
95  {
96  if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
97  if (nbrhs == 0) x = NULL;
98  z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm);
99  }
100 
101  // Convert the matrix to Fortran-style Numbering
102  template <typename MatrixType>
103  void EigenToFortranNumbering (MatrixType& mat)
104  {
105  if ( !(mat.outerIndexPtr()[0]) )
106  {
107  int i;
108  for(i = 0; i <= mat.rows(); ++i)
109  ++mat.outerIndexPtr()[i];
110  for(i = 0; i < mat.nonZeros(); ++i)
111  ++mat.innerIndexPtr()[i];
112  }
113  }
114 
115  // Convert to C-style Numbering
116  template <typename MatrixType>
117  void EigenToCNumbering (MatrixType& mat)
118  {
119  // Check the Numbering
120  if ( mat.outerIndexPtr()[0] == 1 )
121  { // Convert to C-style numbering
122  int i;
123  for(i = 0; i <= mat.rows(); ++i)
124  --mat.outerIndexPtr()[i];
125  for(i = 0; i < mat.nonZeros(); ++i)
126  --mat.innerIndexPtr()[i];
127  }
128  }
129 
130  // Symmetrize the graph of the input matrix
131  // In : The Input matrix to symmetrize the pattern
132  // Out : The output matrix
133  // StrMatTrans : The structural pattern of the transpose of In; It is
134  // used to optimize the future symmetrization with the same matrix pattern
135  // WARNING It is assumed here that successive calls to this routine are done
136  // with matrices having the same pattern.
137  template <typename MatrixType>
138  void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans, bool& hasTranspose)
139  {
140  eigen_assert(In.cols()==In.rows() && " Can only symmetrize the graph of a square matrix");
141  if (!hasTranspose)
142  { //First call to this routine, need to compute the structural pattern of In^T
143  StrMatTrans = In.transpose();
144  // Set the elements of the matrix to zero
145  for (int i = 0; i < StrMatTrans.rows(); i++)
146  {
147  for (typename MatrixType::InnerIterator it(StrMatTrans, i); it; ++it)
148  it.valueRef() = 0.0;
149  }
150  hasTranspose = true;
151  }
152  Out = (StrMatTrans + In).eval();
153  }
154 
155 }
156 
157 // This is the base class to interface with PaStiX functions.
158 // Users should not used this class directly.
159 template <class Derived>
161 {
162  public:
163  typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
165  typedef typename MatrixType::Scalar Scalar;
166  typedef typename MatrixType::RealScalar RealScalar;
167  typedef typename MatrixType::Index Index;
169 
170  public:
171 
173  {
174  m_pastixdata = 0;
175  m_hasTranspose = false;
176  PastixInit();
177  }
178 
180  {
181  PastixDestroy();
182  }
183 
184  // Initialize the Pastix data structure, check the matrix
185  void PastixInit();
186 
187  // Compute the ordering and the symbolic factorization
188  Derived& analyzePattern (MatrixType& mat);
189 
190  // Compute the numerical factorization
191  Derived& factorize (MatrixType& mat);
192 
197  template<typename Rhs>
198  inline const internal::solve_retval<PastixBase, Rhs>
199  solve(const MatrixBase<Rhs>& b) const
200  {
201  eigen_assert(m_isInitialized && "Pastix solver is not initialized.");
202  eigen_assert(rows()==b.rows()
203  && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
204  return internal::solve_retval<PastixBase, Rhs>(*this, b.derived());
205  }
206 
207  template<typename Rhs,typename Dest>
208  bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
209 
211  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
213  {
214  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
215  eigen_assert(rows()==b.rows());
216 
217  // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
218  static const int NbColsAtOnce = 1;
219  int rhsCols = b.cols();
220  int size = b.rows();
222  for(int k=0; k<rhsCols; k+=NbColsAtOnce)
223  {
224  int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
225  tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
226  tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
227  dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
228  }
229  }
230 
231  Derived& derived()
232  {
233  return *static_cast<Derived*>(this);
234  }
235  const Derived& derived() const
236  {
237  return *static_cast<const Derived*>(this);
238  }
239 
246  {
247  return m_iparm;
248  }
249 
254  int& iparm(int idxparam)
255  {
256  return m_iparm(idxparam);
257  }
258 
264  {
265  return m_dparm;
266  }
267 
268 
273  double& dparm(int idxparam)
274  {
275  return m_dparm(idxparam);
276  }
277 
278  inline Index cols() const { return m_size; }
279  inline Index rows() const { return m_size; }
280 
290  {
291  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
292  return m_info;
293  }
294 
299  template<typename Rhs>
300  inline const internal::sparse_solve_retval<PastixBase, Rhs>
301  solve(const SparseMatrixBase<Rhs>& b) const
302  {
303  eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized.");
304  eigen_assert(rows()==b.rows()
305  && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
306  return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived());
307  }
308 
309  protected:
310  // Free all the data allocated by Pastix
312  {
313  eigen_assert(m_initisOk && "The Pastix structure should be allocated first");
314  m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
315  m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
318  }
319 
320  Derived& compute (MatrixType& mat);
321 
327  mutable pastix_data_t *m_pastixdata; // Data structure for pastix
328  mutable SparseMatrix<Scalar, ColMajor> m_mat_null; // An input null matrix
329  mutable Matrix<Scalar, Dynamic,1> m_vec_null; // An input null vector
330  mutable SparseMatrix<Scalar, ColMajor> m_StrMatTrans; // The transpose pattern of the input matrix
331  mutable bool m_hasTranspose; // The transpose of the current matrix has already been computed
332  mutable int m_comm; // The MPI communicator identifier
333  mutable Matrix<Index,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
334  mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
335  mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector
336  mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector
337  mutable int m_ordering; // ordering method to use
338  mutable int m_amalgamation; // level of amalgamation
339  mutable int m_size; // Size of the matrix
340 
341  private:
342  PastixBase(PastixBase& ) {}
343 
344 };
345 
350 template <class Derived>
352 {
353  m_size = 0;
354  m_iparm.resize(IPARM_SIZE);
355  m_dparm.resize(DPARM_SIZE);
356 
357  m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
358  if(m_pastixdata)
359  { // This trick is used to reset the Pastix internal data between successive
360  // calls with (structural) different matrices
361  PastixDestroy();
362  m_pastixdata = 0;
363  m_iparm(IPARM_MODIFY_PARAMETER) = API_YES;
364  m_hasTranspose = false;
365  }
366 
367  m_iparm(IPARM_START_TASK) = API_TASK_INIT;
368  m_iparm(IPARM_END_TASK) = API_TASK_INIT;
369  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
370  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
371  m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data());
372 
373  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
374 
375  // Check the returned error
376  if(m_iparm(IPARM_ERROR_NUMBER)) {
377  m_info = InvalidInput;
378  m_initisOk = false;
379  }
380  else {
381  m_info = Success;
382  m_initisOk = true;
383  }
384 }
385 
386 template <class Derived>
388 {
389  eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
390  typedef typename MatrixType::Scalar Scalar;
391 
392  // Save the size of the current matrix
393  m_size = mat.rows();
394  // Convert the matrix in fortran-style numbering
396  analyzePattern(mat);
397  factorize(mat);
398  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
399  if (m_factorizationIsOk) m_isInitialized = true;
400 
401  //Convert back the matrix -- Is it really necessary here
403 
404  return derived();
405 }
406 
407 
408 template <class Derived>
410 {
411  eigen_assert(m_initisOk && "PastixInit should be called first to set the default parameters");
412  m_size = mat.rows();
413  m_perm.resize(m_size);
414  m_invp.resize(m_size);
415 
416  // Convert the matrix in fortran-style numbering
418 
419  m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
420  m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
421 
422  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
423  mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data());
424 
425  // Check the returned error
426  if(m_iparm(IPARM_ERROR_NUMBER)) {
427  m_info = NumericalIssue;
428  m_analysisIsOk = false;
429  }
430  else {
431  m_info = Success;
432  m_analysisIsOk = true;
433  }
434  return derived();
435 }
436 
437 template <class Derived>
439 {
440  eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
441  m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
442  m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
443  m_size = mat.rows();
444 
445  // Convert the matrix in fortran-style numbering
447 
448  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
449  mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data());
450 
451  // Check the returned error
452  if(m_iparm(IPARM_ERROR_NUMBER)) {
453  m_info = NumericalIssue;
454  m_factorizationIsOk = false;
455  m_isInitialized = false;
456  }
457  else {
458  m_info = Success;
459  m_factorizationIsOk = true;
460  m_isInitialized = true;
461  }
462  return derived();
463 }
464 
465 /* Solve the system */
466 template<typename Base>
467 template<typename Rhs,typename Dest>
469 {
470  eigen_assert(m_isInitialized && "The matrix should be factorized first");
471  EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
472  THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
473  int rhs = 1;
474 
475  x = b; /* on return, x is overwritten by the computed solution */
476 
477  for (int i = 0; i < b.cols(); i++){
478  m_iparm(IPARM_START_TASK) = API_TASK_SOLVE;
479  m_iparm(IPARM_END_TASK) = API_TASK_REFINE;
480  m_iparm(IPARM_RHS_MAKING) = API_RHS_B;
481  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
482  m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
483  }
484  // Check the returned error
485  if(m_iparm(IPARM_ERROR_NUMBER)) {
486  m_info = NumericalIssue;
487  return false;
488  }
489  else {
490  return true;
491  }
492 }
493 
513 template<typename _MatrixType, bool IsStrSym>
514 class PastixLU : public PastixBase< PastixLU<_MatrixType> >
515 {
516  public:
519  typedef typename MatrixType::Scalar Scalar;
521 
522  public:
523  PastixLU():Base() {}
524 
525  PastixLU(const MatrixType& matrix):Base()
526  {
527  compute(matrix);
528  }
534  void compute (const MatrixType& matrix)
535  {
536  // Pastix supports only column-major matrices with a symmetric pattern
537  Base::PastixInit();
538  PaStiXType temp(matrix.rows(), matrix.cols());
539  // Symmetrize the graph of the matrix
540  if (IsStrSym)
541  temp = matrix;
542  else
543  {
544  internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans, m_hasTranspose);
545  }
546  m_iparm[IPARM_SYM] = API_SYM_NO;
547  m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
548  Base::compute(temp);
549  }
555  void analyzePattern(const MatrixType& matrix)
556  {
557 
558  Base::PastixInit();
559  /* Pastix supports only column-major matrices with symmetrized patterns */
560  SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
561  // Symmetrize the graph of the matrix
562  if (IsStrSym)
563  temp = matrix;
564  else
565  {
566  internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
567  }
568 
569  m_iparm(IPARM_SYM) = API_SYM_NO;
570  m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
571  Base::analyzePattern(temp);
572  }
573 
579  void factorize(const MatrixType& matrix)
580  {
581  /* Pastix supports only column-major matrices with symmetrized patterns */
582  SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
583  // Symmetrize the graph of the matrix
584  if (IsStrSym)
585  temp = matrix;
586  else
587  {
588  internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
589  }
590  m_iparm(IPARM_SYM) = API_SYM_NO;
591  m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
592  Base::factorize(temp);
593  }
594  protected:
595  using Base::m_iparm;
596  using Base::m_dparm;
597  using Base::m_StrMatTrans;
598  using Base::m_hasTranspose;
599 
600  private:
601  PastixLU(PastixLU& ) {}
602 };
603 
618 template<typename _MatrixType, int _UpLo>
619 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
620 {
621  public:
624  typedef typename MatrixType::Scalar Scalar;
625  typedef typename MatrixType::Index Index;
626 
627  public:
628  enum { UpLo = _UpLo };
630 
631  PastixLLT(const MatrixType& matrix):Base()
632  {
633  compute(matrix);
634  }
635 
639  void compute (const MatrixType& matrix)
640  {
641  // Pastix supports only lower, column-major matrices
642  Base::PastixInit(); // This is necessary to let PaStiX initialize its data structure between successive calls to compute
643  SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
645  temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
646  m_iparm(IPARM_SYM) = API_SYM_YES;
647  m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
648  Base::compute(temp);
649  }
650 
655  void analyzePattern(const MatrixType& matrix)
656  {
657  Base::PastixInit();
658  // Pastix supports only lower, column-major matrices
659  SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
661  temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
662  m_iparm(IPARM_SYM) = API_SYM_YES;
663  m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
664  Base::analyzePattern(temp);
665  }
669  void factorize(const MatrixType& matrix)
670  {
671  // Pastix supports only lower, column-major matrices
672  SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
674  temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
675  m_iparm(IPARM_SYM) = API_SYM_YES;
676  m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
677  Base::factorize(temp);
678  }
679  protected:
680  using Base::m_iparm;
681 
682  private:
683  PastixLLT(PastixLLT& ) {}
684 };
685 
700 template<typename _MatrixType, int _UpLo>
701 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
702 {
703 public:
706  typedef typename MatrixType::Scalar Scalar;
707  typedef typename MatrixType::Index Index;
708 
709  public:
710  enum { UpLo = _UpLo };
712 
713  PastixLDLT(const MatrixType& matrix):Base()
714  {
715  compute(matrix);
716  }
717 
721  void compute (const MatrixType& matrix)
722  {
724  // Pastix supports only lower, column-major matrices
725  SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
727  temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
728  m_iparm(IPARM_SYM) = API_SYM_YES;
729  m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
730  Base::compute(temp);
731  }
732 
737  void analyzePattern(const MatrixType& matrix)
738  {
740  // Pastix supports only lower, column-major matrices
741  SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
743  temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
744 
745  m_iparm(IPARM_SYM) = API_SYM_YES;
746  m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
747  Base::analyzePattern(temp);
748  }
752  void factorize(const MatrixType& matrix)
753  {
754  // Pastix supports only lower, column-major matrices
755  SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
757  temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
758 
759  m_iparm(IPARM_SYM) = API_SYM_YES;
760  m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
761  Base::factorize(temp);
762  }
763 
764  protected:
765  using Base::m_iparm;
766 
767  private:
768  PastixLDLT(PastixLDLT& ) {}
769 };
770 
771 namespace internal {
772 
773 template<typename _MatrixType, typename Rhs>
774 struct solve_retval<PastixBase<_MatrixType>, Rhs>
775  : solve_retval_base<PastixBase<_MatrixType>, Rhs>
776 {
777  typedef PastixBase<_MatrixType> Dec;
778  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
779 
780  template<typename Dest> void evalTo(Dest& dst) const
781  {
782  dec()._solve(rhs(),dst);
783  }
784 };
785 
786 template<typename _MatrixType, typename Rhs>
787 struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs>
788  : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs>
789 {
790  typedef PastixBase<_MatrixType> Dec;
792 
793  template<typename Dest> void evalTo(Dest& dst) const
794  {
795  dec()._solve_sparse(rhs(),dst);
796  }
797 };
798 
799 } // end namespace internal
800 
801 } // end namespace Eigen
802 
803 #endif