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  {
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 
134  {
135  pardisoRelease();
136  }
137 
138  inline Index cols() const { return m_size; }
139  inline Index rows() const { return m_size; }
140 
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:
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  bool symmetric = std::abs(m_type) < 10;
196  m_iparm[0] = 1; // No solver default
197  m_iparm[1] = 2; // use Metis for the ordering
198  m_iparm[2] = 0; // Reserved. Set to zero. (??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, b is left unchanged
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[36] = 0; // CSR
223  m_iparm[59] = 0; // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core
224 
225  memset(m_pt, 0, sizeof(m_pt));
226  }
227 
228  protected:
229  // cached data to reduce reallocation, etc.
230 
232  {
233  switch(error)
234  {
235  case 0:
236  m_info = Success;
237  break;
238  case -4:
239  case -7:
240  m_info = NumericalIssue;
241  break;
242  default:
243  m_info = InvalidInput;
244  }
245  }
246 
247  mutable SparseMatrixType m_matrix;
249  bool m_analysisIsOk, m_factorizationIsOk;
250  StorageIndex m_type, m_msglvl;
251  mutable void *m_pt[64];
252  mutable ParameterType m_iparm;
253  mutable IntColVectorType m_perm;
255 
256 };
257 
258 template<class Derived>
260 {
261  m_size = a.rows();
262  eigen_assert(a.rows() == a.cols());
263 
264  pardisoRelease();
265  m_perm.setZero(m_size);
266  derived().getMatrix(a);
267 
268  Index error;
269  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size),
270  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
271  m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
272  manageErrorCode(error);
273  m_analysisIsOk = true;
274  m_factorizationIsOk = true;
275  m_isInitialized = true;
276  return derived();
277 }
278 
279 template<class Derived>
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, internal::convert_index<StorageIndex>(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>
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, internal::convert_index<StorageIndex>(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>
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, internal::convert_index<StorageIndex>(m_size),
358  m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
359  m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl,
360  rhs_ptr, x.derived().data());
361 
362  manageErrorCode(error);
363 }
364 
365 
383 template<typename MatrixType>
384 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
385 {
386  protected:
388  typedef typename Base::Scalar Scalar;
389  typedef typename Base::RealScalar RealScalar;
390  using Base::pardisoInit;
391  using Base::m_matrix;
392  friend class PardisoImpl< PardisoLU<MatrixType> >;
393 
394  public:
395 
396  using Base::compute;
397  using Base::solve;
398 
400  : Base()
401  {
402  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
403  }
404 
405  explicit PardisoLU(const MatrixType& matrix)
406  : Base()
407  {
408  pardisoInit(Base::ScalarIsComplex ? 13 : 11);
409  compute(matrix);
410  }
411  protected:
413  {
414  m_matrix = matrix;
415  m_matrix.makeCompressed();
416  }
417 };
418 
438 template<typename MatrixType, int _UpLo>
439 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
440 {
441  protected:
443  typedef typename Base::Scalar Scalar;
444  typedef typename Base::RealScalar RealScalar;
445  using Base::pardisoInit;
446  using Base::m_matrix;
447  friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
448 
449  public:
450 
452  enum { UpLo = _UpLo };
453  using Base::compute;
454 
456  : Base()
457  {
458  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
459  }
460 
461  explicit PardisoLLT(const MatrixType& matrix)
462  : Base()
463  {
464  pardisoInit(Base::ScalarIsComplex ? 4 : 2);
465  compute(matrix);
466  }
467 
468  protected:
469 
471  {
472  // PARDISO supports only upper, row-major matrices
474  m_matrix.resize(matrix.rows(), matrix.cols());
475  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
476  m_matrix.makeCompressed();
477  }
478 };
479 
501 template<typename MatrixType, int Options>
502 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
503 {
504  protected:
506  typedef typename Base::Scalar Scalar;
507  typedef typename Base::RealScalar RealScalar;
508  using Base::pardisoInit;
509  using Base::m_matrix;
510  friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
511 
512  public:
513 
515  using Base::compute;
516  enum { UpLo = Options&(Upper|Lower) };
517 
519  : Base()
520  {
521  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
522  }
523 
524  explicit PardisoLDLT(const MatrixType& matrix)
525  : Base()
526  {
527  pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
528  compute(matrix);
529  }
530 
532  {
533  // PARDISO supports only upper, row-major matrices
535  m_matrix.resize(matrix.rows(), matrix.cols());
536  m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
537  m_matrix.makeCompressed();
538  }
539 };
540 
541 } // end namespace Eigen
542 
543 #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:33
Matrix< StorageIndex, 1, MatrixType::ColsAtCompileTime > IntRowVectorType
Scalar * b
Definition: benchVecAdd.cpp:17
return int(ret)+1
SparseMatrix< Scalar, RowMajor, StorageIndex > SparseMatrixType
Derived & compute(const MatrixType &matrix)
Derived & factorize(const MatrixType &matrix)
Index rows() const
void _solve_impl(const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
PardisoLDLT(const MatrixType &matrix)
A base class for sparse solvers.
int n
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
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:150
Array< StorageIndex, 64, 1, DontAlign > ParameterType
static const Point3 pt(1.0, 2.0, 3.0)
void getMatrix(const MatrixType &matrix)
const unsigned int RowMajorBit
Definition: Constants.h:61
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
Array33i a
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)
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:33
#define eigen_assert(x)
Definition: Macros.h:579
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:34
Base::Scalar Scalar
void manageErrorCode(Index error) const
internal::pardiso_traits< Derived > Traits
ComputationInfo info() const
Reports whether previous computation was successful.
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:39
const int Dynamic
Definition: Constants.h:21
ParameterType & pardisoParameterArray()
void getMatrix(const MatrixType &matrix)
Base::Scalar Scalar
Index cols() const
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:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
IntColVectorType m_perm
SparseMatrixType m_matrix
Definition: pytypes.h:897
PardisoLLT(const MatrixType &matrix)
PardisoLU(const MatrixType &matrix)
Base::StorageIndex StorageIndex


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:43:22