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 IndexType>
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;
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;
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;
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 
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 
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 
135  {
136  pardisoRelease();
137  }
138 
139  inline Index cols() const { return m_size; }
140  inline Index rows() const { return m_size; }
141 
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:
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 
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;
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;
256 
257 };
258 
259 template<class Derived>
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>
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>
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>
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());
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:
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 
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:
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:
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;
454  enum { UpLo = _UpLo };
455  using Base::compute;
456 
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 
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:
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;
517  using Base::compute;
518  enum { UpLo = Options&(Upper|Lower) };
519 
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 
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
StorageIndex m_type
static IndexType run(_MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x)
SCALAR Scalar
Definition: bench_gemm.cpp:46
Matrix< StorageIndex, 1, MatrixType::ColsAtCompileTime > IntRowVectorType
Scalar * b
Definition: benchVecAdd.cpp:17
SparseMatrix< Scalar, RowMajor, StorageIndex > SparseMatrixType
Derived & compute(const MatrixType &matrix)
Derived & factorize(const MatrixType &matrix)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
PardisoLDLT(const MatrixType &matrix)
A base class for sparse solvers.
int n
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixXf MatrixType
Matrix< StorageIndex, MatrixType::RowsAtCompileTime, 1 > IntColVectorType
Traits::StorageIndex StorageIndex
PardisoImpl< PardisoLLT< MatrixType, _UpLo > > Base
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Array< StorageIndex, 64, 1, DontAlign > ParameterType
static const Point3 pt(1.0, 2.0, 3.0)
void getMatrix(const MatrixType &matrix)
void _solve_impl(const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
const unsigned int RowMajorBit
Definition: Constants.h:66
static IndexType run(_MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a, IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x)
Traits::RealScalar RealScalar
Base::RealScalar RealScalar
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Base::StorageIndex StorageIndex
Base::Scalar Scalar
void pardisoInit(int type)
void manageErrorCode(Index error) const
Derived & analyzePattern(const MatrixType &matrix)
ComputationInfo m_info
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
#define eigen_assert(x)
Definition: Macros.h:1037
Traits::Scalar Scalar
PardisoImpl< PardisoLU > Base
idx_t idx_t idx_t idx_t idx_t * perm
Base::RealScalar RealScalar
A sparse direct LU factorization and solver based on the PARDISO library.
#define NULL
Definition: ccolamd.c:609
Matrix< Scalar, Dynamic, 1 > VectorType
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Base::Scalar Scalar
ComputationInfo info() const
Reports whether previous computation was successful.
internal::pardiso_traits< Derived > Traits
void getMatrix(const MatrixType &matrix)
Base::RealScalar RealScalar
SparseSolverBase< Derived > Base
PardisoImpl< PardisoLDLT< MatrixType, Options > > Base
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
static double error
Definition: testRot3.cpp:37
const int Dynamic
Definition: Constants.h:22
ParameterType & pardisoParameterArray()
void getMatrix(const MatrixType &matrix)
Base::Scalar Scalar
ParameterType m_iparm
Traits::MatrixType MatrixType
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
void resize(Index newSize)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define abs(x)
Definition: datatypes.h:17
ComputationInfo
Definition: Constants.h:440
Index cols() const
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
IntColVectorType m_perm
Index rows() const
SparseMatrixType m_matrix
Definition: pytypes.h:1370
PardisoLLT(const MatrixType &matrix)
PardisoLU(const MatrixType &matrix)
Base::StorageIndex StorageIndex


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:11