Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.99 (git rev 199c5f2b47eb1f8e5a2d20e60f07e97cd95a6ba6)
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:
102  typedef SparseSolverBase<Derived> Base;
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;
114  typedef SparseMatrix<Scalar,RowMajor,StorageIndex> SparseMatrixType;
115  typedef Matrix<Scalar,Dynamic,1> VectorType;
116  typedef Matrix<StorageIndex, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
117  typedef Matrix<StorageIndex, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
118  typedef Array<StorageIndex,64,1,DontAlign> ParameterType;
119  enum {
120  ScalarIsComplex = NumTraits<Scalar>::IsComplex,
121  ColsAtCompileTime = Dynamic,
122  MaxColsAtCompileTime = Dynamic
123  };
124 
125  PardisoImpl()
126  : m_analysisIsOk(false), m_factorizationIsOk(false)
127  {
128  eigen_assert((sizeof(StorageIndex) >= sizeof(_INTEGER_t) && sizeof(StorageIndex) <= 8) && "Non-supported index type");
129  m_iparm.setZero();
130  m_msglvl = 0; // No output
131  m_isInitialized = false;
132  }
133 
134  ~PardisoImpl()
135  {
136  pardisoRelease();
137  }
138 
139  inline Index cols() const { return m_size; }
140  inline Index rows() const { return m_size; }
141 
147  ComputationInfo info() const
148  {
149  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
150  return m_info;
151  }
152 
156  ParameterType& pardisoParameterArray()
157  {
158  return m_iparm;
159  }
160 
167  Derived& analyzePattern(const MatrixType& matrix);
168 
175  Derived& factorize(const MatrixType& matrix);
176 
177  Derived& compute(const MatrixType& matrix);
178 
179  template<typename Rhs,typename Dest>
180  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
181 
182  protected:
183  void pardisoRelease()
184  {
185  if(m_isInitialized) // Factorization ran at least once
186  {
187  internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, -1, internal::convert_index<StorageIndex>(m_size),0, 0, 0, m_perm.data(), 0,
188  m_iparm.data(), m_msglvl, NULL, NULL);
189  m_isInitialized = false;
190  }
191  }
192 
193  void pardisoInit(int type)
194  {
195  m_type = type;
196  bool symmetric = std::abs(m_type) < 10;
197  m_iparm[0] = 1; // No solver default
198  m_iparm[1] = 2; // use Metis for the ordering
199  m_iparm[2] = 0; // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??)
200  m_iparm[3] = 0; // No iterative-direct algorithm
201  m_iparm[4] = 0; // No user fill-in reducing permutation
202  m_iparm[5] = 0; // Write solution into x, b is left unchanged
203  m_iparm[6] = 0; // Not in use
204  m_iparm[7] = 2; // Max numbers of iterative refinement steps
205  m_iparm[8] = 0; // Not in use
206  m_iparm[9] = 13; // Perturb the pivot elements with 1E-13
207  m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
208  m_iparm[11] = 0; // Not in use
209  m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric).
210  // Try m_iparm[12] = 1 in case of inappropriate accuracy
211  m_iparm[13] = 0; // Output: Number of perturbed pivots
212  m_iparm[14] = 0; // Not in use
213  m_iparm[15] = 0; // Not in use
214  m_iparm[16] = 0; // Not in use
215  m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
216  m_iparm[18] = -1; // Output: Mflops for LU factorization
217  m_iparm[19] = 0; // Output: Numbers of CG Iterations
218 
219  m_iparm[20] = 0; // 1x1 pivoting
220  m_iparm[26] = 0; // No matrix checker
221  m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
222  m_iparm[34] = 1; // C indexing
223  m_iparm[36] = 0; // CSR
224  m_iparm[59] = 0; // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core
225 
226  memset(m_pt, 0, sizeof(m_pt));
227  }
228 
229  protected:
230  // cached data to reduce reallocation, etc.
231 
232  void manageErrorCode(Index error) const
233  {
234  switch(error)
235  {
236  case 0:
237  m_info = Success;
238  break;
239  case -4:
240  case -7:
241  m_info = NumericalIssue;
242  break;
243  default:
244  m_info = InvalidInput;
245  }
246  }
247 
248  mutable SparseMatrixType m_matrix;
249  mutable ComputationInfo m_info;
250  bool m_analysisIsOk, m_factorizationIsOk;
251  StorageIndex m_type, m_msglvl;
252  mutable void *m_pt[64];
253  mutable ParameterType m_iparm;
254  mutable IntColVectorType m_perm;
255  Index m_size;
256 
257 };
258 
259 template<class Derived>
260 Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
261 {
262  m_size = a.rows();
263  eigen_assert(a.rows() == a.cols());
264 
265  pardisoRelease();
266  m_perm.setZero(m_size);
267  derived().getMatrix(a);
268 
269  Index error;
270  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size),
271  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
272  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
273  manageErrorCode(error);
274  m_analysisIsOk = true;
275  m_factorizationIsOk = true;
276  m_isInitialized = true;
277  return derived();
278 }
279 
280 template<class Derived>
281 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
282 {
283  m_size = a.rows();
284  eigen_assert(m_size == a.cols());
285 
286  pardisoRelease();
287  m_perm.setZero(m_size);
288  derived().getMatrix(a);
289 
290  Index error;
291  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, internal::convert_index<StorageIndex>(m_size),
292  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
293  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
294 
295  manageErrorCode(error);
296  m_analysisIsOk = true;
297  m_factorizationIsOk = false;
298  m_isInitialized = true;
299  return derived();
300 }
301 
302 template<class Derived>
303 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
304 {
305  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
306  eigen_assert(m_size == a.rows() && m_size == a.cols());
307 
308  derived().getMatrix(a);
309 
310  Index error;
311  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, internal::convert_index<StorageIndex>(m_size),
312  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
313  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
314 
315  manageErrorCode(error);
316  m_factorizationIsOk = true;
317  return derived();
318 }
319 
320 template<class Derived>
321 template<typename BDerived,typename XDerived>
322 void PardisoImpl<Derived>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
323 {
324  if(m_iparm[0] == 0) // Factorization was not computed
325  {
326  m_info = InvalidInput;
327  return;
328  }
329 
330  //Index n = m_matrix.rows();
331  Index nrhs = Index(b.cols());
332  eigen_assert(m_size==b.rows());
333  eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
334  eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
335  eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
336 
337 
338 // switch (transposed) {
339 // case SvNoTrans : m_iparm[11] = 0 ; break;
340 // case SvTranspose : m_iparm[11] = 2 ; break;
341 // case SvAdjoint : m_iparm[11] = 1 ; break;
342 // default:
343 // //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n";
344 // m_iparm[11] = 0;
345 // }
346 
347  Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
348  Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
349 
350  // Pardiso cannot solve in-place
351  if(rhs_ptr == x.derived().data())
352  {
353  tmp = b;
354  rhs_ptr = tmp.data();
355  }
356 
357  Index error;
358  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, internal::convert_index<StorageIndex>(m_size),
359  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
360  m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl,
361  rhs_ptr, x.derived().data());
362 
363  manageErrorCode(error);
364 }
365 
366 
384 template<typename MatrixType>
385 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
386 {
387  protected:
388  typedef PardisoImpl<PardisoLU> Base;
389  using Base::pardisoInit;
390  using Base::m_matrix;
391  friend class PardisoImpl< PardisoLU<MatrixType> >;
392 
393  public:
394 
395  typedef typename Base::Scalar Scalar;
396  typedef typename Base::RealScalar RealScalar;
397 
398  using Base::compute;
399  using Base::solve;
400 
401  PardisoLU()
402  : Base()
403  {
404  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
405  }
406 
407  explicit PardisoLU(const MatrixType& matrix)
408  : Base()
409  {
410  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
411  compute(matrix);
412  }
413  protected:
414  void getMatrix(const MatrixType& matrix)
415  {
416  m_matrix = matrix;
417  m_matrix.makeCompressed();
418  }
419 };
420 
440 template<typename MatrixType, int _UpLo>
441 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
442 {
443  protected:
444  typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
445  using Base::pardisoInit;
446  using Base::m_matrix;
447  friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
448 
449  public:
450 
451  typedef typename Base::Scalar Scalar;
452  typedef typename Base::RealScalar RealScalar;
453  typedef typename Base::StorageIndex StorageIndex;
454  enum { UpLo = _UpLo };
455  using Base::compute;
456 
457  PardisoLLT()
458  : Base()
459  {
460  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
461  }
462 
463  explicit PardisoLLT(const MatrixType& matrix)
464  : Base()
465  {
466  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
467  compute(matrix);
468  }
469 
470  protected:
471 
472  void getMatrix(const MatrixType& matrix)
473  {
474  // PARDISO supports only upper, row-major matrices
476  m_matrix.resize(matrix.rows(), matrix.cols());
477  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
478  m_matrix.makeCompressed();
479  }
480 };
481 
503 template<typename MatrixType, int Options>
504 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
505 {
506  protected:
507  typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
508  using Base::pardisoInit;
509  using Base::m_matrix;
510  friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
511 
512  public:
513 
514  typedef typename Base::Scalar Scalar;
515  typedef typename Base::RealScalar RealScalar;
516  typedef typename Base::StorageIndex StorageIndex;
517  using Base::compute;
518  enum { UpLo = Options&(Upper|Lower) };
519 
520  PardisoLDLT()
521  : Base()
522  {
523  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
524  }
525 
526  explicit PardisoLDLT(const MatrixType& matrix)
527  : Base()
528  {
529  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
530  compute(matrix);
531  }
532 
533  void getMatrix(const MatrixType& matrix)
534  {
535  // PARDISO supports only upper, row-major matrices
537  m_matrix.resize(matrix.rows(), matrix.cols());
538  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
539  m_matrix.makeCompressed();
540  }
541 };
542 
543 } // end namespace Eigen
544 
545 #endif // EIGEN_PARDISOSUPPORT_H
@ Flags
Definition: DenseBase.h:165
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:505
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:442
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:386
void resize(Index newSize)
Definition: PermutationMatrix.h:125
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
ComputationInfo
Definition: Constants.h:440
@ Symmetric
Definition: Constants.h:227
@ Lower
Definition: Constants.h:209
@ Upper
Definition: Constants.h:211
@ NumericalIssue
Definition: Constants.h:444
@ InvalidInput
Definition: Constants.h:449
@ Success
Definition: Constants.h:442
const unsigned int RowMajorBit
Definition: Constants.h:66
Namespace containing all symbols from the Eigen library.
Definition: Core:134
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22