Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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
35namespace Eigen {
36
37template<typename _MatrixType> class PardisoLU;
38template<typename _MatrixType, int Options=Upper> class PardisoLLT;
39template<typename _MatrixType, int Options=Upper> class PardisoLDLT;
40
41namespace 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
98template<class Derived>
99class 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 : 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
134 ~PardisoImpl()
135 {
136 pardisoRelease();
137 }
138
139 inline Index cols() const { return m_size; }
140 inline Index rows() const { return m_size; }
141
147 ComputationInfo info() const
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:
183 void pardisoRelease()
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
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
259template<class Derived>
260Derived& 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
280template<class Derived>
281Derived& 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
302template<class Derived>
303Derived& 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
320template<class Derived>
321template<typename BDerived,typename XDerived>
322void 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
384template<typename MatrixType>
385class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
386{
387 protected:
388 typedef PardisoImpl<PardisoLU> Base;
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
401 PardisoLU()
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:
414 void getMatrix(const MatrixType& matrix)
415 {
416 m_matrix = matrix;
417 m_matrix.makeCompressed();
418 }
419};
420
440template<typename MatrixType, int _UpLo>
441class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
442{
443 protected:
444 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
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;
453 typedef typename Base::StorageIndex StorageIndex;
454 enum { UpLo = _UpLo };
455 using Base::compute;
456
457 PardisoLLT()
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
472 void getMatrix(const MatrixType& matrix)
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
503template<typename MatrixType, int Options>
504class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
505{
506 protected:
507 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
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;
516 typedef typename Base::StorageIndex StorageIndex;
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
533 void getMatrix(const MatrixType& matrix)
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
A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:505
A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:442
A sparse direct LU factorization and solver based on the PARDISO library.
Definition: PardisoSupport.h:386
void resize(Index newSize)
Definition: PermutationMatrix.h:125
Permutation matrix.
Definition: PermutationMatrix.h:298
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SparseSolverBase.h:88
ComputationInfo
Definition: Constants.h:440
@ Symmetric
Definition: Constants.h:227
@ Lower
Definition: Constants.h:209
@ Upper
Definition: Constants.h:211
@ NumericalIssue
Definition: Constants.h:444
@ InvalidInput
Definition: Constants.h:449
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22