PardisoSupport.h
Go to the documentation of this file.
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>
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;
111  enum {
113  };
114 
116  {
117  eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type");
118  m_iparm.setZero();
119  m_msglvl = 0; // No output
120  m_initialized = false;
121  }
122 
124  {
125  pardisoRelease();
126  }
127 
128  inline Index cols() const { return m_size; }
129  inline Index rows() const { return m_size; }
130 
137  {
138  eigen_assert(m_initialized && "Decomposition is not initialized.");
139  return m_info;
140  }
141 
146  {
147  return m_iparm;
148  }
149 
156  Derived& analyzePattern(const MatrixType& matrix);
157 
164  Derived& factorize(const MatrixType& matrix);
165 
166  Derived& compute(const MatrixType& matrix);
167 
172  template<typename Rhs>
173  inline const internal::solve_retval<PardisoImpl, Rhs>
174  solve(const MatrixBase<Rhs>& b) const
175  {
176  eigen_assert(m_initialized && "Pardiso solver is not initialized.");
177  eigen_assert(rows()==b.rows()
178  && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
179  return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived());
180  }
181 
186  template<typename Rhs>
187  inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
188  solve(const SparseMatrixBase<Rhs>& b) const
189  {
190  eigen_assert(m_initialized && "Pardiso solver is not initialized.");
191  eigen_assert(rows()==b.rows()
192  && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
193  return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived());
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 
205  template<typename BDerived, typename XDerived>
206  bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const;
207 
209  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
211  {
212  eigen_assert(m_size==b.rows());
213 
214  // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
215  static const int NbColsAtOnce = 4;
216  int rhsCols = b.cols();
217  int size = b.rows();
218  // Pardiso cannot solve in-place,
219  // so we need two temporaries
222  for(int k=0; k<rhsCols; k+=NbColsAtOnce)
223  {
224  int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
225  tmp_rhs.leftCols(actualCols) = b.middleCols(k,actualCols);
226  tmp_res.leftCols(actualCols) = derived().solve(tmp_rhs.leftCols(actualCols));
227  dest.middleCols(k,actualCols) = tmp_res.leftCols(actualCols).sparseView();
228  }
229  }
230 
231  protected:
233  {
234  if(m_initialized) // Factorization ran at least once
235  {
236  internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
237  m_iparm.data(), m_msglvl, 0, 0);
238  }
239  }
240 
241  void pardisoInit(int type)
242  {
243  m_type = type;
244  bool symmetric = abs(m_type) < 10;
245  m_iparm[0] = 1; // No solver default
246  m_iparm[1] = 3; // use Metis for the ordering
247  m_iparm[2] = 1; // Numbers of processors, value of OMP_NUM_THREADS
248  m_iparm[3] = 0; // No iterative-direct algorithm
249  m_iparm[4] = 0; // No user fill-in reducing permutation
250  m_iparm[5] = 0; // Write solution into x
251  m_iparm[6] = 0; // Not in use
252  m_iparm[7] = 2; // Max numbers of iterative refinement steps
253  m_iparm[8] = 0; // Not in use
254  m_iparm[9] = 13; // Perturb the pivot elements with 1E-13
255  m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
256  m_iparm[11] = 0; // Not in use
257  m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric).
258  // Try m_iparm[12] = 1 in case of inappropriate accuracy
259  m_iparm[13] = 0; // Output: Number of perturbed pivots
260  m_iparm[14] = 0; // Not in use
261  m_iparm[15] = 0; // Not in use
262  m_iparm[16] = 0; // Not in use
263  m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
264  m_iparm[18] = -1; // Output: Mflops for LU factorization
265  m_iparm[19] = 0; // Output: Numbers of CG Iterations
266 
267  m_iparm[20] = 0; // 1x1 pivoting
268  m_iparm[26] = 0; // No matrix checker
269  m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
270  m_iparm[34] = 1; // C indexing
271  m_iparm[59] = 1; // Automatic switch between In-Core and Out-of-Core modes
272  }
273 
274  protected:
275  // cached data to reduce reallocation, etc.
276 
277  void manageErrorCode(Index error)
278  {
279  switch(error)
280  {
281  case 0:
282  m_info = Success;
283  break;
284  case -4:
285  case -7:
287  break;
288  default:
290  }
291  }
292 
297  mutable void *m_pt[64];
301 
302  private:
304 };
305 
306 template<class Derived>
308 {
309  m_size = a.rows();
310  eigen_assert(a.rows() == a.cols());
311 
312  pardisoRelease();
313  memset(m_pt, 0, sizeof(m_pt));
314  m_perm.setZero(m_size);
315  derived().getMatrix(a);
316 
317  Index error;
318  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
319  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
320  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
321 
322  manageErrorCode(error);
323  m_analysisIsOk = true;
324  m_factorizationIsOk = true;
325  m_initialized = true;
326  return derived();
327 }
328 
329 template<class Derived>
331 {
332  m_size = a.rows();
333  eigen_assert(m_size == a.cols());
334 
335  pardisoRelease();
336  memset(m_pt, 0, sizeof(m_pt));
337  m_perm.setZero(m_size);
338  derived().getMatrix(a);
339 
340  Index error;
341  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
342  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
343  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
344 
345  manageErrorCode(error);
346  m_analysisIsOk = true;
347  m_factorizationIsOk = false;
348  m_initialized = true;
349  return derived();
350 }
351 
352 template<class Derived>
354 {
355  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
356  eigen_assert(m_size == a.rows() && m_size == a.cols());
357 
358  derived().getMatrix(a);
359 
360  Index error;
361  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
362  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
363  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
364 
365  manageErrorCode(error);
366  m_factorizationIsOk = true;
367  return derived();
368 }
369 
370 template<class Base>
371 template<typename BDerived,typename XDerived>
373 {
374  if(m_iparm[0] == 0) // Factorization was not computed
375  return false;
376 
377  //Index n = m_matrix.rows();
378  Index nrhs = Index(b.cols());
379  eigen_assert(m_size==b.rows());
380  eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
381  eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
382  eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
383 
384 
385 // switch (transposed) {
386 // case SvNoTrans : m_iparm[11] = 0 ; break;
387 // case SvTranspose : m_iparm[11] = 2 ; break;
388 // case SvAdjoint : m_iparm[11] = 1 ; break;
389 // default:
390 // //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n";
391 // m_iparm[11] = 0;
392 // }
393 
394  Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
396 
397  // Pardiso cannot solve in-place
398  if(rhs_ptr == x.derived().data())
399  {
400  tmp = b;
401  rhs_ptr = tmp.data();
402  }
403 
404  Index error;
405  error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
406  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
407  m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
408  rhs_ptr, x.derived().data());
409 
410  return error==0;
411 }
412 
413 
426 template<typename MatrixType>
427 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
428 {
429  protected:
431  typedef typename Base::Scalar Scalar;
432  typedef typename Base::RealScalar RealScalar;
433  using Base::pardisoInit;
434  using Base::m_matrix;
435  friend class PardisoImpl< PardisoLU<MatrixType> >;
436 
437  public:
438 
439  using Base::compute;
440  using Base::solve;
441 
443  : Base()
444  {
446  }
447 
448  PardisoLU(const MatrixType& matrix)
449  : Base()
450  {
452  compute(matrix);
453  }
454  protected:
455  void getMatrix(const MatrixType& matrix)
456  {
457  m_matrix = matrix;
458  }
459 
460  private:
461  PardisoLU(PardisoLU& ) {}
462 };
463 
478 template<typename MatrixType, int _UpLo>
479 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
480 {
481  protected:
483  typedef typename Base::Scalar Scalar;
484  typedef typename Base::Index Index;
485  typedef typename Base::RealScalar RealScalar;
486  using Base::pardisoInit;
487  using Base::m_matrix;
488  friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
489 
490  public:
491 
492  enum { UpLo = _UpLo };
493  using Base::compute;
494  using Base::solve;
495 
497  : Base()
498  {
500  }
501 
502  PardisoLLT(const MatrixType& matrix)
503  : Base()
504  {
506  compute(matrix);
507  }
508 
509  protected:
510 
511  void getMatrix(const MatrixType& matrix)
512  {
513  // PARDISO supports only upper, row-major matrices
515  m_matrix.resize(matrix.rows(), matrix.cols());
516  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
517  }
518 
519  private:
520  PardisoLLT(PardisoLLT& ) {}
521 };
522 
539 template<typename MatrixType, int Options>
540 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
541 {
542  protected:
544  typedef typename Base::Scalar Scalar;
545  typedef typename Base::Index Index;
546  typedef typename Base::RealScalar RealScalar;
547  using Base::pardisoInit;
548  using Base::m_matrix;
549  friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
550 
551  public:
552 
553  using Base::compute;
554  using Base::solve;
555  enum { UpLo = Options&(Upper|Lower) };
556 
558  : Base()
559  {
560  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
561  }
562 
563  PardisoLDLT(const MatrixType& matrix)
564  : Base()
565  {
566  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
567  compute(matrix);
568  }
569 
570  void getMatrix(const MatrixType& matrix)
571  {
572  // PARDISO supports only upper, row-major matrices
574  m_matrix.resize(matrix.rows(), matrix.cols());
575  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
576  }
577 
578  private:
580 };
581 
582 namespace internal {
583 
584 template<typename _Derived, typename Rhs>
585 struct solve_retval<PardisoImpl<_Derived>, Rhs>
586  : solve_retval_base<PardisoImpl<_Derived>, Rhs>
587 {
588  typedef PardisoImpl<_Derived> Dec;
589  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
590 
591  template<typename Dest> void evalTo(Dest& dst) const
592  {
593  dec()._solve(rhs(),dst);
594  }
595 };
596 
597 template<typename Derived, typename Rhs>
598 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
599  : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
600 {
601  typedef PardisoImpl<Derived> Dec;
603 
604  template<typename Dest> void evalTo(Dest& dst) const
605  {
606  dec().derived()._solve_sparse(rhs(),dst);
607  }
608 };
609 
610 } // end namespace internal
611 
612 } // end namespace Eigen
613 
614 #endif // EIGEN_PARDISOSUPPORT_H