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