Eigen  3.2.92
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 IndexType>
44  struct pardiso_run_selector
45  {
46  static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a,
47  IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x)
48  {
49  IndexType 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 IndexType;
58  static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a,
59  IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x)
60  {
61  IndexType 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::StorageIndex StorageIndex;
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::StorageIndex StorageIndex;
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::StorageIndex StorageIndex;
94  };
95 
96 } // end namespace internal
97 
98 template<class Derived>
99 class PardisoImpl : public SparseSolverBase<Derived>
100 {
101  protected:
103  using Base::derived;
104  using Base::m_isInitialized;
105 
106  typedef internal::pardiso_traits<Derived> Traits;
107  public:
108  using Base::_solve_impl;
109 
110  typedef typename Traits::MatrixType MatrixType;
111  typedef typename Traits::Scalar Scalar;
112  typedef typename Traits::RealScalar RealScalar;
113  typedef typename Traits::StorageIndex StorageIndex;
119  enum {
120  ScalarIsComplex = NumTraits<Scalar>::IsComplex,
121  ColsAtCompileTime = Dynamic,
122  MaxColsAtCompileTime = Dynamic
123  };
124 
125  PardisoImpl()
126  {
127  eigen_assert((sizeof(StorageIndex) >= sizeof(_INTEGER_t) && sizeof(StorageIndex) <= 8) && "Non-supported index type");
128  m_iparm.setZero();
129  m_msglvl = 0; // No output
130  m_isInitialized = false;
131  }
132 
133  ~PardisoImpl()
134  {
135  pardisoRelease();
136  }
137 
138  inline Index cols() const { return m_size; }
139  inline Index rows() const { return m_size; }
140 
146  ComputationInfo info() const
147  {
148  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
149  return m_info;
150  }
151 
155  ParameterType& pardisoParameterArray()
156  {
157  return m_iparm;
158  }
159 
166  Derived& analyzePattern(const MatrixType& matrix);
167 
174  Derived& factorize(const MatrixType& matrix);
175 
176  Derived& compute(const MatrixType& matrix);
177 
178  template<typename Rhs,typename Dest>
179  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
180 
181  protected:
182  void pardisoRelease()
183  {
184  if(m_isInitialized) // Factorization ran at least once
185  {
186  internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, -1, m_size,0, 0, 0, m_perm.data(), 0,
187  m_iparm.data(), m_msglvl, NULL, NULL);
188  m_isInitialized = false;
189  }
190  }
191 
192  void pardisoInit(int type)
193  {
194  m_type = type;
195  bool symmetric = std::abs(m_type) < 10;
196  m_iparm[0] = 1; // No solver default
197  m_iparm[1] = 3; // use Metis for the ordering
198  m_iparm[2] = 1; // Numbers of processors, value of OMP_NUM_THREADS
199  m_iparm[3] = 0; // No iterative-direct algorithm
200  m_iparm[4] = 0; // No user fill-in reducing permutation
201  m_iparm[5] = 0; // Write solution into x
202  m_iparm[6] = 0; // Not in use
203  m_iparm[7] = 2; // Max numbers of iterative refinement steps
204  m_iparm[8] = 0; // Not in use
205  m_iparm[9] = 13; // Perturb the pivot elements with 1E-13
206  m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
207  m_iparm[11] = 0; // Not in use
208  m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric).
209  // Try m_iparm[12] = 1 in case of inappropriate accuracy
210  m_iparm[13] = 0; // Output: Number of perturbed pivots
211  m_iparm[14] = 0; // Not in use
212  m_iparm[15] = 0; // Not in use
213  m_iparm[16] = 0; // Not in use
214  m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
215  m_iparm[18] = -1; // Output: Mflops for LU factorization
216  m_iparm[19] = 0; // Output: Numbers of CG Iterations
217 
218  m_iparm[20] = 0; // 1x1 pivoting
219  m_iparm[26] = 0; // No matrix checker
220  m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
221  m_iparm[34] = 1; // C indexing
222  m_iparm[59] = 1; // Automatic switch between In-Core and Out-of-Core modes
223 
224  memset(m_pt, 0, sizeof(m_pt));
225  }
226 
227  protected:
228  // cached data to reduce reallocation, etc.
229 
230  void manageErrorCode(Index error) const
231  {
232  switch(error)
233  {
234  case 0:
235  m_info = Success;
236  break;
237  case -4:
238  case -7:
239  m_info = NumericalIssue;
240  break;
241  default:
242  m_info = InvalidInput;
243  }
244  }
245 
246  mutable SparseMatrixType m_matrix;
247  mutable ComputationInfo m_info;
248  bool m_analysisIsOk, m_factorizationIsOk;
249  Index m_type, m_msglvl;
250  mutable void *m_pt[64];
251  mutable ParameterType m_iparm;
252  mutable IntColVectorType m_perm;
253  Index m_size;
254 
255 };
256 
257 template<class Derived>
258 Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
259 {
260  m_size = a.rows();
261  eigen_assert(a.rows() == a.cols());
262 
263  pardisoRelease();
264  m_perm.setZero(m_size);
265  derived().getMatrix(a);
266 
267  Index error;
268  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, m_size,
269  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
270  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
271 
272  manageErrorCode(error);
273  m_analysisIsOk = true;
274  m_factorizationIsOk = true;
275  m_isInitialized = true;
276  return derived();
277 }
278 
279 template<class Derived>
280 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
281 {
282  m_size = a.rows();
283  eigen_assert(m_size == a.cols());
284 
285  pardisoRelease();
286  m_perm.setZero(m_size);
287  derived().getMatrix(a);
288 
289  Index error;
290  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, m_size,
291  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
292  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
293 
294  manageErrorCode(error);
295  m_analysisIsOk = true;
296  m_factorizationIsOk = false;
297  m_isInitialized = true;
298  return derived();
299 }
300 
301 template<class Derived>
302 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
303 {
304  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
305  eigen_assert(m_size == a.rows() && m_size == a.cols());
306 
307  derived().getMatrix(a);
308 
309  Index error;
310  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, m_size,
311  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
312  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
313 
314  manageErrorCode(error);
315  m_factorizationIsOk = true;
316  return derived();
317 }
318 
319 template<class Derived>
320 template<typename BDerived,typename XDerived>
321 void PardisoImpl<Derived>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
322 {
323  if(m_iparm[0] == 0) // Factorization was not computed
324  {
325  m_info = InvalidInput;
326  return;
327  }
328 
329  //Index n = m_matrix.rows();
330  Index nrhs = Index(b.cols());
331  eigen_assert(m_size==b.rows());
332  eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
333  eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
334  eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
335 
336 
337 // switch (transposed) {
338 // case SvNoTrans : m_iparm[11] = 0 ; break;
339 // case SvTranspose : m_iparm[11] = 2 ; break;
340 // case SvAdjoint : m_iparm[11] = 1 ; break;
341 // default:
342 // //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n";
343 // m_iparm[11] = 0;
344 // }
345 
346  Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
348 
349  // Pardiso cannot solve in-place
350  if(rhs_ptr == x.derived().data())
351  {
352  tmp = b;
353  rhs_ptr = tmp.data();
354  }
355 
356  Index error;
357  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, m_size,
358  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
359  m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
360  rhs_ptr, x.derived().data());
361 
362  manageErrorCode(error);
363 }
364 
365 
380 template<typename MatrixType>
381 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
382 {
383  protected:
384  typedef PardisoImpl<PardisoLU> Base;
385  typedef typename Base::Scalar Scalar;
386  typedef typename Base::RealScalar RealScalar;
387  using Base::pardisoInit;
388  using Base::m_matrix;
389  friend class PardisoImpl< PardisoLU<MatrixType> >;
390 
391  public:
392 
393  using Base::compute;
394  using Base::solve;
395 
396  PardisoLU()
397  : Base()
398  {
399  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
400  }
401 
402  explicit PardisoLU(const MatrixType& matrix)
403  : Base()
404  {
405  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
406  compute(matrix);
407  }
408  protected:
409  void getMatrix(const MatrixType& matrix)
410  {
411  m_matrix = matrix;
412  m_matrix.makeCompressed();
413  }
414 };
415 
432 template<typename MatrixType, int _UpLo>
433 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
434 {
435  protected:
436  typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
437  typedef typename Base::Scalar Scalar;
438  typedef typename Base::RealScalar RealScalar;
439  using Base::pardisoInit;
440  using Base::m_matrix;
441  friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
442 
443  public:
444 
445  typedef typename Base::StorageIndex StorageIndex;
446  enum { UpLo = _UpLo };
447  using Base::compute;
448 
449  PardisoLLT()
450  : Base()
451  {
452  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
453  }
454 
455  explicit PardisoLLT(const MatrixType& matrix)
456  : Base()
457  {
458  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
459  compute(matrix);
460  }
461 
462  protected:
463 
464  void getMatrix(const MatrixType& matrix)
465  {
466  // PARDISO supports only upper, row-major matrices
468  m_matrix.resize(matrix.rows(), matrix.cols());
469  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
470  m_matrix.makeCompressed();
471  }
472 };
473 
492 template<typename MatrixType, int Options>
493 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
494 {
495  protected:
496  typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
497  typedef typename Base::Scalar Scalar;
498  typedef typename Base::RealScalar RealScalar;
499  using Base::pardisoInit;
500  using Base::m_matrix;
501  friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
502 
503  public:
504 
505  typedef typename Base::StorageIndex StorageIndex;
506  using Base::compute;
507  enum { UpLo = Options&(Upper|Lower) };
508 
509  PardisoLDLT()
510  : Base()
511  {
512  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
513  }
514 
515  explicit PardisoLDLT(const MatrixType& matrix)
516  : Base()
517  {
518  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
519  compute(matrix);
520  }
521 
522  void getMatrix(const MatrixType& matrix)
523  {
524  // PARDISO supports only upper, row-major matrices
526  m_matrix.resize(matrix.rows(), matrix.cols());
527  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
528  m_matrix.makeCompressed();
529  }
530 };
531 
532 } // end namespace Eigen
533 
534 #endif // EIGEN_PARDISOSUPPORT_H
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Definition: LDLT.h:16
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
const unsigned int RowMajorBit
Definition: Constants.h:61
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:38
Definition: Constants.h:204
Definition: Constants.h:434
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
Definition: Constants.h:439
Definition: Constants.h:432
Definition: Eigen_Colamd.h:54
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:39
void resize(Index newSize)
Definition: PermutationMatrix.h:137
Definition: Constants.h:222
Definition: Constants.h:206
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48