Eigen  3.2.10
PardisoSupport.h
1 /*
2  Copyright (c) 2011, Intel Corporation. All rights reserved.
3 
4  Redistribution and use in source and binary forms, with or without modification,
5  are permitted provided that the following conditions are met:
6 
7  * Redistributions of source code must retain the above copyright notice, this
8  list of conditions and the following disclaimer.
9  * Redistributions in binary form must reproduce the above copyright notice,
10  this list of conditions and the following disclaimer in the documentation
11  and/or other materials provided with the distribution.
12  * Neither the name of Intel Corporation nor the names of its contributors may
13  be used to endorse or promote products derived from this software without
14  specific prior written permission.
15 
16  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
20  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
23  ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 
27  ********************************************************************************
28  * Content : Eigen bindings to Intel(R) MKL PARDISO
29  ********************************************************************************
30 */
31 
32 #ifndef EIGEN_PARDISOSUPPORT_H
33 #define EIGEN_PARDISOSUPPORT_H
34 
35 namespace Eigen {
36 
37 template<typename _MatrixType> class PardisoLU;
38 template<typename _MatrixType, int Options=Upper> class PardisoLLT;
39 template<typename _MatrixType, int Options=Upper> class PardisoLDLT;
40 
41 namespace internal
42 {
43  template<typename Index>
44  struct pardiso_run_selector
45  {
46  static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a,
47  Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
48  {
49  Index error = 0;
50  ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
51  return error;
52  }
53  };
54  template<>
55  struct pardiso_run_selector<long long int>
56  {
57  typedef long long int Index;
58  static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a,
59  Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
60  {
61  Index error = 0;
62  ::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
63  return error;
64  }
65  };
66 
67  template<class Pardiso> struct pardiso_traits;
68 
69  template<typename _MatrixType>
70  struct pardiso_traits< PardisoLU<_MatrixType> >
71  {
72  typedef _MatrixType MatrixType;
73  typedef typename _MatrixType::Scalar Scalar;
74  typedef typename _MatrixType::RealScalar RealScalar;
75  typedef typename _MatrixType::Index Index;
76  };
77 
78  template<typename _MatrixType, int Options>
79  struct pardiso_traits< PardisoLLT<_MatrixType, Options> >
80  {
81  typedef _MatrixType MatrixType;
82  typedef typename _MatrixType::Scalar Scalar;
83  typedef typename _MatrixType::RealScalar RealScalar;
84  typedef typename _MatrixType::Index Index;
85  };
86 
87  template<typename _MatrixType, int Options>
88  struct pardiso_traits< PardisoLDLT<_MatrixType, Options> >
89  {
90  typedef _MatrixType MatrixType;
91  typedef typename _MatrixType::Scalar Scalar;
92  typedef typename _MatrixType::RealScalar RealScalar;
93  typedef typename _MatrixType::Index Index;
94  };
95 
96 }
97 
98 template<class Derived>
99 class PardisoImpl
100 {
101  typedef internal::pardiso_traits<Derived> Traits;
102  public:
103  typedef typename Traits::MatrixType MatrixType;
104  typedef typename Traits::Scalar Scalar;
105  typedef typename Traits::RealScalar RealScalar;
106  typedef typename Traits::Index Index;
112  enum {
113  ScalarIsComplex = NumTraits<Scalar>::IsComplex
114  };
115 
116  PardisoImpl()
117  {
118  eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type");
119  m_iparm.setZero();
120  m_msglvl = 0; // No output
121  m_initialized = false;
122  }
123 
124  ~PardisoImpl()
125  {
126  pardisoRelease();
127  }
128 
129  inline Index cols() const { return m_size; }
130  inline Index rows() const { return m_size; }
131 
137  ComputationInfo info() const
138  {
139  eigen_assert(m_initialized && "Decomposition is not initialized.");
140  return m_info;
141  }
142 
146  ParameterType& pardisoParameterArray()
147  {
148  return m_iparm;
149  }
150 
157  Derived& analyzePattern(const MatrixType& matrix);
158 
165  Derived& factorize(const MatrixType& matrix);
166 
167  Derived& compute(const MatrixType& matrix);
168 
173  template<typename Rhs>
174  inline const internal::solve_retval<PardisoImpl, Rhs>
175  solve(const MatrixBase<Rhs>& b) const
176  {
177  eigen_assert(m_initialized && "Pardiso solver is not initialized.");
178  eigen_assert(rows()==b.rows()
179  && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
180  return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived());
181  }
182 
187  template<typename Rhs>
188  inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
189  solve(const SparseMatrixBase<Rhs>& b) const
190  {
191  eigen_assert(m_initialized && "Pardiso solver is not initialized.");
192  eigen_assert(rows()==b.rows()
193  && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
194  return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived());
195  }
196 
197  Derived& derived()
198  {
199  return *static_cast<Derived*>(this);
200  }
201  const Derived& derived() const
202  {
203  return *static_cast<const Derived*>(this);
204  }
205 
206  template<typename BDerived, typename XDerived>
207  bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const;
208 
209  protected:
210  void pardisoRelease()
211  {
212  if(m_initialized) // Factorization ran at least once
213  {
214  internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
215  m_iparm.data(), m_msglvl, 0, 0);
216  }
217  }
218 
219  void pardisoInit(int type)
220  {
221  m_type = type;
222  bool symmetric = std::abs(m_type) < 10;
223  m_iparm[0] = 1; // No solver default
224  m_iparm[1] = 2; // use Metis for the ordering
225  m_iparm[2] = 0; // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??)
226  m_iparm[3] = 0; // No iterative-direct algorithm
227  m_iparm[4] = 0; // No user fill-in reducing permutation
228  m_iparm[5] = 0; // Write solution into x, b is left unchanged
229  m_iparm[6] = 0; // Not in use
230  m_iparm[7] = 2; // Max numbers of iterative refinement steps
231  m_iparm[8] = 0; // Not in use
232  m_iparm[9] = 13; // Perturb the pivot elements with 1E-13
233  m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
234  m_iparm[11] = 0; // Not in use
235  m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric).
236  // Try m_iparm[12] = 1 in case of inappropriate accuracy
237  m_iparm[13] = 0; // Output: Number of perturbed pivots
238  m_iparm[14] = 0; // Not in use
239  m_iparm[15] = 0; // Not in use
240  m_iparm[16] = 0; // Not in use
241  m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
242  m_iparm[18] = -1; // Output: Mflops for LU factorization
243  m_iparm[19] = 0; // Output: Numbers of CG Iterations
244 
245  m_iparm[20] = 0; // 1x1 pivoting
246  m_iparm[26] = 0; // No matrix checker
247  m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
248  m_iparm[34] = 1; // C indexing
249  m_iparm[36] = 0; // CSR
250  m_iparm[59] = 0; // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core
251 
252  memset(m_pt, 0, sizeof(m_pt));
253  }
254 
255  protected:
256  // cached data to reduce reallocation, etc.
257 
258  void manageErrorCode(Index error)
259  {
260  switch(error)
261  {
262  case 0:
263  m_info = Success;
264  break;
265  case -4:
266  case -7:
267  m_info = NumericalIssue;
268  break;
269  default:
270  m_info = InvalidInput;
271  }
272  }
273 
274  mutable SparseMatrixType m_matrix;
275  ComputationInfo m_info;
276  bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
277  Index m_type, m_msglvl;
278  mutable void *m_pt[64];
279  mutable ParameterType m_iparm;
280  mutable IntColVectorType m_perm;
281  Index m_size;
282 
283  private:
284  PardisoImpl(PardisoImpl &) {}
285 };
286 
287 template<class Derived>
288 Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
289 {
290  m_size = a.rows();
291  eigen_assert(a.rows() == a.cols());
292 
293  pardisoRelease();
294  memset(m_pt, 0, sizeof(m_pt));
295  m_perm.setZero(m_size);
296  derived().getMatrix(a);
297 
298  Index error;
299  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
300  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
301  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
302 
303  manageErrorCode(error);
304  m_analysisIsOk = true;
305  m_factorizationIsOk = true;
306  m_initialized = true;
307  return derived();
308 }
309 
310 template<class Derived>
311 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
312 {
313  m_size = a.rows();
314  eigen_assert(m_size == a.cols());
315 
316  pardisoRelease();
317  memset(m_pt, 0, sizeof(m_pt));
318  m_perm.setZero(m_size);
319  derived().getMatrix(a);
320 
321  Index error;
322  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
323  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
324  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
325 
326  manageErrorCode(error);
327  m_analysisIsOk = true;
328  m_factorizationIsOk = false;
329  m_initialized = true;
330  return derived();
331 }
332 
333 template<class Derived>
334 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
335 {
336  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
337  eigen_assert(m_size == a.rows() && m_size == a.cols());
338 
339  derived().getMatrix(a);
340 
341  Index error;
342  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
343  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
344  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
345 
346  manageErrorCode(error);
347  m_factorizationIsOk = true;
348  return derived();
349 }
350 
351 template<class Base>
352 template<typename BDerived,typename XDerived>
353 bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
354 {
355  if(m_iparm[0] == 0) // Factorization was not computed
356  return false;
357 
358  //Index n = m_matrix.rows();
359  Index nrhs = Index(b.cols());
360  eigen_assert(m_size==b.rows());
361  eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
362  eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
363  eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
364 
365 
366 // switch (transposed) {
367 // case SvNoTrans : m_iparm[11] = 0 ; break;
368 // case SvTranspose : m_iparm[11] = 2 ; break;
369 // case SvAdjoint : m_iparm[11] = 1 ; break;
370 // default:
371 // //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n";
372 // m_iparm[11] = 0;
373 // }
374 
375  Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
377 
378  // Pardiso cannot solve in-place
379  if(rhs_ptr == x.derived().data())
380  {
381  tmp = b;
382  rhs_ptr = tmp.data();
383  }
384 
385  Index error;
386  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
387  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
388  m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
389  rhs_ptr, x.derived().data());
390  return error==0;
391 }
392 
393 
409 template<typename MatrixType>
410 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
411 {
412  protected:
413  typedef PardisoImpl< PardisoLU<MatrixType> > Base;
414  typedef typename Base::Scalar Scalar;
415  typedef typename Base::RealScalar RealScalar;
416  using Base::pardisoInit;
417  using Base::m_matrix;
418  friend class PardisoImpl< PardisoLU<MatrixType> >;
419 
420  public:
421 
422  using Base::compute;
423  using Base::solve;
424 
425  PardisoLU()
426  : Base()
427  {
428  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
429  }
430 
431  PardisoLU(const MatrixType& matrix)
432  : Base()
433  {
434  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
435  compute(matrix);
436  }
437  protected:
438  void getMatrix(const MatrixType& matrix)
439  {
440  m_matrix = matrix;
441  }
442 
443  private:
444  PardisoLU(PardisoLU& ) {}
445 };
446 
464 template<typename MatrixType, int _UpLo>
465 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
466 {
467  protected:
468  typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
469  typedef typename Base::Scalar Scalar;
470  typedef typename Base::Index Index;
471  typedef typename Base::RealScalar RealScalar;
472  using Base::pardisoInit;
473  using Base::m_matrix;
474  friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
475 
476  public:
477 
478  enum { UpLo = _UpLo };
479  using Base::compute;
480  using Base::solve;
481 
482  PardisoLLT()
483  : Base()
484  {
485  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
486  }
487 
488  PardisoLLT(const MatrixType& matrix)
489  : Base()
490  {
491  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
492  compute(matrix);
493  }
494 
495  protected:
496 
497  void getMatrix(const MatrixType& matrix)
498  {
499  // PARDISO supports only upper, row-major matrices
501  m_matrix.resize(matrix.rows(), matrix.cols());
502  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
503  }
504 
505  private:
506  PardisoLLT(PardisoLLT& ) {}
507 };
508 
528 template<typename MatrixType, int Options>
529 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
530 {
531  protected:
532  typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
533  typedef typename Base::Scalar Scalar;
534  typedef typename Base::Index Index;
535  typedef typename Base::RealScalar RealScalar;
536  using Base::pardisoInit;
537  using Base::m_matrix;
538  friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
539 
540  public:
541 
542  using Base::compute;
543  using Base::solve;
544  enum { UpLo = Options&(Upper|Lower) };
545 
546  PardisoLDLT()
547  : Base()
548  {
549  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
550  }
551 
552  PardisoLDLT(const MatrixType& matrix)
553  : Base()
554  {
555  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
556  compute(matrix);
557  }
558 
559  void getMatrix(const MatrixType& matrix)
560  {
561  // PARDISO supports only upper, row-major matrices
563  m_matrix.resize(matrix.rows(), matrix.cols());
564  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
565  }
566 
567  private:
568  PardisoLDLT(PardisoLDLT& ) {}
569 };
570 
571 namespace internal {
572 
573 template<typename _Derived, typename Rhs>
574 struct solve_retval<PardisoImpl<_Derived>, Rhs>
575  : solve_retval_base<PardisoImpl<_Derived>, Rhs>
576 {
577  typedef PardisoImpl<_Derived> Dec;
578  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
579 
580  template<typename Dest> void evalTo(Dest& dst) const
581  {
582  dec()._solve(rhs(),dst);
583  }
584 };
585 
586 template<typename Derived, typename Rhs>
587 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
588  : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
589 {
590  typedef PardisoImpl<Derived> Dec;
591  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
592 
593  template<typename Dest> void evalTo(Dest& dst) const
594  {
595  this->defaultEvalTo(dst);
596  }
597 };
598 
599 } // end namespace internal
600 
601 } // end namespace Eigen
602 
603 #endif // EIGEN_PARDISOSUPPORT_H
Definition: Constants.h:378
Definition: LDLT.h:16
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Definition: Constants.h:185
Index rows() const
Definition: SparseMatrixBase.h:160
Definition: Constants.h:167
void resize(Index newSize)
Definition: PermutationMatrix.h:142
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:38
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:239
Derived & derived()
Definition: EigenBase.h:34
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
Definition: Constants.h:383
Definition: Eigen_Colamd.h:50
Definition: Constants.h:169
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:39
Definition: Constants.h:376
const unsigned int RowMajorBit
Definition: Constants.h:53
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48