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