Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.3.8
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  {
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, internal::convert_index<StorageIndex>(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  EIGEN_USING_STD(abs);
196  bool symmetric = 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  typedef typename Base::Scalar Scalar;
390  typedef typename Base::RealScalar RealScalar;
391  using Base::pardisoInit;
392  using Base::m_matrix;
393  friend class PardisoImpl< PardisoLU<MatrixType> >;
394 
395  public:
396 
397  using Base::compute;
398  using Base::solve;
399 
400  PardisoLU()
401  : Base()
402  {
403  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
404  }
405 
406  explicit PardisoLU(const MatrixType& matrix)
407  : Base()
408  {
409  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
410  compute(matrix);
411  }
412  protected:
413  void getMatrix(const MatrixType& matrix)
414  {
415  m_matrix = matrix;
416  m_matrix.makeCompressed();
417  }
418 };
419 
439 template<typename MatrixType, int _UpLo>
440 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
441 {
442  protected:
443  typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
444  typedef typename Base::Scalar Scalar;
445  typedef typename Base::RealScalar RealScalar;
446  using Base::pardisoInit;
447  using Base::m_matrix;
448  friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
449 
450  public:
451 
452  typedef typename Base::StorageIndex StorageIndex;
453  enum { UpLo = _UpLo };
454  using Base::compute;
455 
456  PardisoLLT()
457  : Base()
458  {
459  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
460  }
461 
462  explicit PardisoLLT(const MatrixType& matrix)
463  : Base()
464  {
465  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
466  compute(matrix);
467  }
468 
469  protected:
470 
471  void getMatrix(const MatrixType& matrix)
472  {
473  // PARDISO supports only upper, row-major matrices
474  PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
475  m_matrix.resize(matrix.rows(), matrix.cols());
476  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
477  m_matrix.makeCompressed();
478  }
479 };
480 
502 template<typename MatrixType, int Options>
503 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
504 {
505  protected:
506  typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
507  typedef typename Base::Scalar Scalar;
508  typedef typename Base::RealScalar RealScalar;
509  using Base::pardisoInit;
510  using Base::m_matrix;
511  friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
512 
513  public:
514 
515  typedef typename Base::StorageIndex StorageIndex;
516  using Base::compute;
517  enum { UpLo = Options&(Upper|Lower) };
518 
519  PardisoLDLT()
520  : Base()
521  {
522  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
523  }
524 
525  explicit PardisoLDLT(const MatrixType& matrix)
526  : Base()
527  {
528  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
529  compute(matrix);
530  }
531 
532  void getMatrix(const MatrixType& matrix)
533  {
534  // PARDISO supports only upper, row-major matrices
535  PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
536  m_matrix.resize(matrix.rows(), matrix.cols());
537  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
538  m_matrix.makeCompressed();
539  }
540 };
541 
542 } // end namespace Eigen
543 
544 #endif // EIGEN_PARDISOSUPPORT_H
Eigen::NumericalIssue
@ NumericalIssue
Definition: Constants.h:434
Eigen::Symmetric
@ Symmetric
Definition: Constants.h:222
Eigen
Namespace containing all symbols from the Eigen library.
Definition: Core:309
Eigen::PardisoLDLT
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:39
Eigen::DenseBase::Flags
@ Flags
Definition: DenseBase.h:160
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:61
Eigen::Upper
@ Upper
Definition: Constants.h:206
Eigen::Success
@ Success
Definition: Constants.h:432
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:21
Eigen::Lower
@ Lower
Definition: Constants.h:204
Eigen::abs
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Eigen::PardisoLU
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:37
Eigen::InvalidInput
@ InvalidInput
Definition: Constants.h:439
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:430
Eigen::PardisoLLT
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:38
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33