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