b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2pardiso_solver.H
1//------------------------------------------------------------------------
2// b2pardiso_solver.H --
3//
4// written by Mathias Doreille
5//
6// (c) 2013 SMR Engineering & Development SA
7// 2502 Bienne, Switzerland
8//
9// All Rights Reserved. Proprietary source code. The contents of
10// this file may not be disclosed to third parties, copied or
11// duplicated in any form, in whole or in part, without the prior
12// written permission of SMR.
13//------------------------------------------------------------------------
14
15#ifndef B2PARDISO_SOLVER_H_
16#define B2PARDISO_SOLVER_H_
17
18#include "/opt/intel/composerxe/mkl/include/mkl_pardiso.h"
19#include "/opt/intel/composerxe/mkl/include/mkl_types.h"
20#include "b2sparse_solver.H"
21#include "utils/b2logging.H"
22
23#define PARDISO64
24
25#ifndef PARDISO64
26#include <vector>
27#endif
28#include <algorithm>
29
30namespace b2000 { namespace b2linalg {
31
32template <typename T>
33class PARDISO_LDLt_seq_sparse_direct_solver : public LDLt_sparse_solver<T> {
34public:
35 PARDISO_LDLt_seq_sparse_direct_solver(bool definit_positive = false) {
36 UnimplementedError() << "The PARDISO library is not enabled." << THROW;
37 }
38
39 void init(
40 size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const T* value,
41 const int connectivity, const Dictionary& dictionary) {}
42
43 void update_value() {}
44
45 void resolve(
46 size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx,
47 char left_or_right = ' ') {
48 UnimplementedError() << THROW;
49 }
50};
51
52template <typename T>
53class PARDISO_LDLt_seq_extension_sparse_direct_solver : public LDLt_extension_sparse_solver<T> {
54public:
55 PARDISO_LDLt_seq_extension_sparse_direct_solver(bool definit_positive = false) {
56 UnimplementedError() << "The PARDISO library is not enabled." << THROW;
57 }
58
59 void init(
60 size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const T* value,
61 size_t s_ext, const int connectivity, const Dictionary& dictionary) {}
62
63 void update_value() {}
64
65 void resolve(
66 size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx, const T* ma = 0,
67 const T* mb = 0, const T* mc = 0, char left_or_right = ' ') {
68 UnimplementedError() << THROW;
69 }
70};
71
72template <typename T>
73class PARDISO_LU_seq_sparse_direct_solver : public LU_sparse_solver<T> {
74public:
75 PARDISO_LU_seq_sparse_direct_solver() {
76 UnimplementedError() << "The PARDISO library is not enabled." << THROW;
77 }
78
79 void init(
80 size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const T* value,
81 const int connectivity, const Dictionary& dictionary) {}
82
83 void update_value() {}
84
85 void resolve(
86 size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx,
87 char left_or_right = ' ') {
88 UnimplementedError() << THROW;
89 }
90};
91
92template <typename T>
93class PARDISO_LU_seq_extension_sparse_direct_solver : public LU_extension_sparse_solver<T> {
94public:
95 PARDISO_LU_seq_extension_sparse_direct_solver() {
96 UnimplementedError() << "The PARDISO library is not enabled." << THROW;
97 }
98
99 void init(
100 size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const T* value,
101 size_t s_ext, const int connectivity, const Dictionary& dictionary) {}
102
103 void update_value() {}
104
105 void resolve(
106 size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx, const T* ma = 0,
107 const T* mb = 0, const T* mc = 0, char left_or_right = ' ') {
108 UnimplementedError() << THROW;
109 }
110};
111
112template <typename T>
113inline void decode_pardiso_param(const Dictionary& d, T* iparam) {
114 std::string icntl = d.get_string("PARDISO_IPARAM", "");
115 if (icntl != "") {
116 std::istringstream iss(icntl);
117 for (;;) {
118 int k, v;
119 char c;
120 if (iss >> k && iss >> c && iss >> v) {
121 if (c != ':' || (k <= 0 && k > 64)) {
122 Exception() << "Cannot parse the PARDISO_IPARAM analysis directive " << icntl
123 << THROW;
124 }
125 iparam[k - 1] = v;
126 } else {
127 Exception() << "Cannot parse the PARDISO_IPARAM analysis directive " << icntl
128 << THROW;
129 }
130 iss >> c;
131 if (!iss) { break; }
132 if (c != ',') {
133 Exception() << "Cannot parse the PARDISO_IPARAM analysis directive " << icntl
134 << THROW;
135 }
136 }
137 }
138}
139
140template <typename T>
141inline std::string decode_pardiso_error_message(T error) {
142 switch (error) {
143 case -1:
144 return "input inconsistent";
145 case -2:
146 return "not enough memory";
147 case -3:
148 return "reordering problem";
149 case -4:
150 return "zero pivot, numerical factorization or iterative refinement problem";
151 case -5:
152 return "unclassified (internal) error";
153 case -6:
154 return "reordering failed (matrix types 11 and 13 only)";
155 case -7:
156 return "diagonal matrix is singular";
157 case -8:
158 return "32-bit integer overflow problem";
159 case -9:
160 return "not enough memory for OOC";
161 case -10:
162 return "problems with opening OOC temporary files";
163 case -11:
164 return "read/write problems with the OOC data file";
165 default:
166 break;
167 }
168 return "Unknown error code";
169}
170
171template <>
172class PARDISO_LDLt_seq_sparse_direct_solver<double> : public LDLt_sparse_solver<double> {
173public:
174 PARDISO_LDLt_seq_sparse_direct_solver(bool definit_positive = false)
175 : updated_value(false), n(0), colptr(0), row(0), avals(0) {
176 std::fill_n(pt, 64, (void*)(0));
177 mtype = definit_positive ? 2 : -2;
178#ifdef PARDISO64
179 int imtype = definit_positive ? 2 : -2;
180 int iiparm[64];
181 pardisoinit(pt, &imtype, iiparm);
182 std::copy(iiparm, iiparm + 64, iparm);
183#else
184 pardisoinit(pt, &mtype, iparm);
185#endif
186 logging::Logger logger =
187 logging::get_logger("linear_algebra.sparse_symmetric_solver.pardiso");
188 msglvl = logger.is_enabled_for(logging::debug) ? 1 : 0;
189 }
190
191 ~PARDISO_LDLt_seq_sparse_direct_solver() {
192#ifdef PARDISO64
193 long long int maxfct = 1;
194 long long int mnum = 1;
195 long long int phase = -1; /* Release internal memory. */
196 long long int idum;
197 long long int nrhs = 1;
198 double ddum; /* Double dummy */
199 long long int error;
200 pardiso_64(
201 pt, &maxfct, &mnum, &mtype, &phase, &n, &ddum, colptr, row, &idum, &nrhs, iparm,
202 &msglvl, &ddum, &ddum, &error);
203#else
204 int maxfct = 1;
205 int mnum = 1;
206 int phase = -1; /* Release internal memory. */
207 int idum;
208 int nrhs = 1;
209 double ddum; /* Double dummy */
210 int error;
211 pardiso(
212 pt, &maxfct, &mnum, &mtype, &phase, &n, &ddum, &colptr[0], &row[0], &idum, &nrhs,
213 iparm, &msglvl, &ddum, &ddum, &error);
214#endif
215 if (error) {
216 Exception() << "The PARDISO solver return error " << error << ", "
217 << decode_pardiso_error_message(error) << THROW;
218 }
219 }
220
221 void init(
222 size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const double* value,
223 const int connectivity, const Dictionary& dictionary) {
224 if (s == 0) { return; }
225 decode_pardiso_param(dictionary, iparm);
226
227 n = s;
228 avals = const_cast<double*>(value);
229#ifdef PARDISO64
230 iparm[34] = 1;
231 colptr = reinterpret_cast<long long int*>(const_cast<size_t*>(colind));
232 row = reinterpret_cast<long long int*>(const_cast<size_t*>(rowind));
233 long long int maxfct = 1;
234 long long int mnum = 1;
235 long long int phase = 12;
236 long long int idum;
237 long long int nrhs = 0;
238 double ddum;
239 long long int error;
240 pardiso_64(
241 pt, &maxfct, &mnum, &mtype, &phase, &n, avals, colptr, row, &idum, &nrhs, iparm,
242 &msglvl, &ddum, &ddum, &error);
243#else
244 colptr.resize(s + 1);
245 row.resize(nnz);
246 {
247 size_t ptr_in = colind[0];
248 size_t ptr_out = 0;
249 colptr[0] = 1;
250 for (size_t j = 1; j <= s; ++j) {
251 colptr[j] = colind[j] + 1;
252 const size_t ptr_end = colind[j];
253 for (; ptr_in != ptr_end; ++ptr_in, ++ptr_out) {
254 row[ptr_out] = rowind[ptr_in] + 1;
255 }
256 }
257 }
258
259 int maxfct = 1;
260 int mnum = 1;
261 int phase = 12;
262 int idum;
263 int nrhs = 0;
264 double ddum;
265 int error;
266 pardiso(
267 pt, &maxfct, &mnum, &mtype, &phase, &n, avals, &colptr[0], &row[0], &idum, &nrhs,
268 iparm, &msglvl, &ddum, &ddum, &error);
269#endif
270 if (error) {
271 Exception() << "The PARDISO solver return error " << error << ", "
272 << decode_pardiso_error_message(error) << THROW;
273 }
274 updated_value = false;
275 }
276
277 void update_value() { updated_value = true; }
278
279 void resolve(
280 size_t s, size_t nrhs, const double* b, size_t ldb, double* x, size_t ldx,
281 char left_or_right = ' ') {
282 if (s == 0) { return; }
283#ifdef PARDISO64
284 long long int maxfct = 1;
285 long long int mnum = 1;
286 long long int phase = updated_value ? 23 : 33;
287 long long int inrhs = nrhs;
288 long long int idum;
289 long long int error;
290 pardiso_64(
291 pt, &maxfct, &mnum, &mtype, &phase, &n, avals, colptr, row, &idum, &inrhs, iparm,
292 &msglvl, const_cast<double*>(b), x, &error);
293#else
294 int maxfct = 1;
295 int mnum = 1;
296 int phase = updated_value ? 23 : 33;
297 int inrhs = nrhs;
298 int idum;
299 int error;
300 pardiso(
301 pt, &maxfct, &mnum, &mtype, &phase, &n, avals, &colptr[0], &row[0], &idum, &inrhs,
302 iparm, &msglvl, const_cast<double*>(b), x, &error);
303#endif
304 if (error) {
305 Exception() << "The PARDISO solver return error " << error << ", "
306 << decode_pardiso_error_message(error) << THROW;
307 }
308 updated_value = false;
309 }
310
311protected:
312 void* pt[64]; /* Pointer to a storage structure needed by PARDISO */
313 bool updated_value;
314#ifdef PARDISO64
315 long long int n;
316 long long int* colptr;
317 long long int* row;
318 long long int mtype;
319 long long int msglvl;
320 long long int iparm[64]; /* integer parameters for PARDISO */
321#else
322 int n;
323 std::vector<int> colptr;
324 std::vector<int> row;
325 int mtype;
326 int msglvl;
327 int iparm[64]; /* integer parameters for PARDISO */
328#endif
329 double* avals;
330};
331
332template <>
333class PARDISO_LDLt_seq_extension_sparse_direct_solver<double>
334 : public LDLt_extension_sparse_solver<double>,
335 public PARDISO_LDLt_seq_sparse_direct_solver<double> {
336public:
337 PARDISO_LDLt_seq_extension_sparse_direct_solver(bool definit_positive = false)
338 : PARDISO_LDLt_seq_sparse_direct_solver<double>(definit_positive), div(0) {}
339
340 void init(
341 size_t size_, size_t nnz_, const size_t* colind_, const size_t* rowind_,
342 const double* value_, size_t size_ext_, const int connectivity,
343 const Dictionary& dictionary) {
344 if (size_ext_ != 1) { UnimplementedError() << THROW; }
345
346 PARDISO_LDLt_seq_sparse_direct_solver<double>::init(
347 size_, nnz_, colind_, rowind_, value_, connectivity, dictionary);
348 m_ab.resize(size_ * 2);
349 }
350
351 void update_value() { PARDISO_LDLt_seq_sparse_direct_solver<double>::update_value(); }
352
353 void resolve(
354 size_t s, size_t nrhs, const double* b, size_t ldb, double* x, size_t ldx,
355 const double* ma_ = 0, const double* mb_ = 0, const double* mc_ = 0,
356 char left_or_right = ' ') {
357 if (mc_ != 0) {
358 std::copy(ma_, ma_ + s - 1, &m_ab[0]);
359 std::copy(mb_, mb_ + s - 1, &m_ab[s - 1]);
360 PARDISO_LDLt_seq_sparse_direct_solver<double>::resolve(
361 s - 1, 2, &m_ab[0], s - 1, &m_ab[0], s - 1);
362 div = 1 / (*mc_ - blas::dot(s - 1, ma_, 1, &m_ab[s - 1], 1));
363 }
364
365 for (size_t i = 0; i != nrhs; ++i) {
366 const double x2 = x[ldx * i + s - 1] =
367 div * (b[ldb * i + s - 1] - blas::dot(s - 1, b + ldb * i, 1, &m_ab[s - 1], 1));
368 PARDISO_LDLt_seq_sparse_direct_solver<double>::resolve(
369 s - 1, 1, b + ldb * i, ldb, x + ldx * i, ldx);
370 blas::axpy(s - 1, -x2, &m_ab[0], 1, x + ldx * i, 1);
371 }
372 }
373
374protected:
375 std::vector<double> m_ab;
376 double div;
377};
378
379template <>
380class PARDISO_LU_seq_sparse_direct_solver<double> : public LU_sparse_solver<double> {
381public:
382 PARDISO_LU_seq_sparse_direct_solver()
383 : updated_value(false), n(0), colptr(0), row(0), avals(0) {
384 std::fill_n(pt, 64, (void*)(0));
385 int mtype = 11;
386#ifdef PARDISO64
387 int iiparm[64];
388 pardisoinit(pt, &mtype, iiparm);
389 std::copy(iiparm, iiparm + 64, iparm);
390#else
391 pardisoinit(pt, &mtype, iparm);
392#endif
393 logging::Logger logger =
394 logging::get_logger("linear_algebra.sparse_symmetric_solver.pardiso");
395 msglvl = logger.is_enabled_for(logging::debug) ? 1 : 0;
396 }
397
398 ~PARDISO_LU_seq_sparse_direct_solver() {
399#ifdef PARDISO64
400 long long int maxfct = 1;
401 long long int mnum = 1;
402 long long int mtype = 11;
403 long long int phase = -1; /* Release internal memory. */
404 long long int idum;
405 long long int nrhs = 1;
406 double ddum; /* Double dummy */
407 long long int error;
408 pardiso_64(
409 pt, &maxfct, &mnum, &mtype, &phase, &n, &ddum, colptr, row, &idum, &nrhs, iparm,
410 &msglvl, &ddum, &ddum, &error);
411#else
412 int maxfct = 1;
413 int mnum = 1;
414 int mtype = 11;
415 int phase = -1; /* Release internal memory. */
416 int idum;
417 int nrhs = 1;
418 double ddum; /* Double dummy */
419 int error;
420 pardiso(
421 pt, &maxfct, &mnum, &mtype, &phase, &n, &ddum, &colptr[0], &row[0], &idum, &nrhs,
422 iparm, &msglvl, &ddum, &ddum, &error);
423#endif
424 if (error) {
425 Exception() << "The PARDISO solver return error " << error << ", "
426 << decode_pardiso_error_message(error) << THROW;
427 }
428 }
429
430 void init(
431 size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const double* value,
432 const int connectivity, const Dictionary& dictionary) {
433 if (s == 0) { return; }
434 decode_pardiso_param(dictionary, iparm);
435
436 n = s;
437 avals = const_cast<double*>(value);
438#ifdef PARDISO64
439 iparm[34] = 1;
440 colptr = reinterpret_cast<long long int*>(const_cast<size_t*>(colind));
441 row = reinterpret_cast<long long int*>(const_cast<size_t*>(rowind));
442 long long int maxfct = 1;
443 long long int mnum = 1;
444 long long int mtype = 11;
445 long long int phase = 12;
446 long long int idum;
447 long long int nrhs = 0;
448 double ddum;
449 long long int error;
450 pardiso_64(
451 pt, &maxfct, &mnum, &mtype, &phase, &n, avals, colptr, row, &idum, &nrhs, iparm,
452 &msglvl, &ddum, &ddum, &error);
453#else
454 colptr.resize(s + 1);
455 row.resize(nnz);
456 {
457 size_t ptr_in = colind[0];
458 size_t ptr_out = 0;
459 colptr[0] = 1;
460 for (size_t j = 1; j <= s; ++j) {
461 colptr[j] = colind[j] + 1;
462 const size_t ptr_end = colind[j];
463 for (; ptr_in != ptr_end; ++ptr_in, ++ptr_out) {
464 row[ptr_out] = rowind[ptr_in] + 1;
465 }
466 }
467 }
468
469 int maxfct = 1;
470 int mnum = 1;
471 int mtype = 11;
472 int phase = 12;
473 int idum;
474 int nrhs = 0;
475 double ddum;
476 int error;
477 pardiso(
478 pt, &maxfct, &mnum, &mtype, &phase, &n, avals, &colptr[0], &row[0], &idum, &nrhs,
479 iparm, &msglvl, &ddum, &ddum, &error);
480#endif
481 if (error) {
482 Exception() << "The PARDISO solver return error " << error << ", "
483 << decode_pardiso_error_message(error) << THROW;
484 }
485 updated_value = false;
486 }
487
488 void update_value() { updated_value = true; }
489
490 void resolve(
491 size_t s, size_t nrhs, const double* b, size_t ldb, double* x, size_t ldx,
492 char left_or_right = ' ') {
493 if (s == 0) { return; }
494#ifdef PARDISO64
495 long long int maxfct = 1;
496 long long int mnum = 1;
497 long long int mtype = 11;
498 long long int phase = updated_value ? 23 : 33;
499 long long int inrhs = nrhs;
500 long long int idum;
501 long long int error;
502 pardiso_64(
503 pt, &maxfct, &mnum, &mtype, &phase, &n, avals, colptr, row, &idum, &inrhs, iparm,
504 &msglvl, const_cast<double*>(b), x, &error);
505#else
506 int maxfct = 1;
507 int mnum = 1;
508 int mtype = 11;
509 int phase = updated_value ? 23 : 33;
510 int inrhs = nrhs;
511 int idum;
512 int error;
513 pardiso(
514 pt, &maxfct, &mnum, &mtype, &phase, &n, avals, &colptr[0], &row[0], &idum, &inrhs,
515 iparm, &msglvl, const_cast<double*>(b), x, &error);
516#endif
517 if (error) {
518 Exception() << "The PARDISO solver return error " << error << ", "
519 << decode_pardiso_error_message(error) << THROW;
520 }
521 updated_value = false;
522 }
523
524protected:
525 void* pt[64]; /* Pointer to a storage structure needed by PARDISO */
526 bool updated_value;
527#ifdef PARDISO64
528 long long int n;
529 long long int* colptr;
530 long long int* row;
531 long long int msglvl;
532 long long int iparm[64]; /* integer parameters for PARDISO */
533#else
534 int n;
535 std::vector<int> colptr;
536 std::vector<int> row;
537 int msglvl;
538 int iparm[64]; /* integer parameters for PARDISO */
539#endif
540 double* avals;
541};
542
543template <>
544class PARDISO_LU_seq_extension_sparse_direct_solver<double>
545 : public LU_extension_sparse_solver<double>,
546 public PARDISO_LU_seq_sparse_direct_solver<double> {
547public:
548 void init(
549 size_t size_, size_t nnz_, const size_t* colind_, const size_t* rowind_,
550 const double* value_, size_t size_ext_, const int connectivity,
551 const Dictionary& dictionary) {
552 if (size_ext_ != 1) { UnimplementedError() << THROW; }
553
554 PARDISO_LU_seq_sparse_direct_solver<double>::init(
555 size_, nnz_, colind_, rowind_, value_, connectivity, dictionary);
556 m_ab.resize(size_ * 2);
557 }
558
559 void update_value() { PARDISO_LU_seq_sparse_direct_solver<double>::update_value(); }
560
561 void resolve(
562 size_t s, size_t nrhs, const double* b, size_t ldb, double* x, size_t ldx,
563 const double* ma_ = 0, const double* mb_ = 0, const double* mc_ = 0,
564 char left_or_right = ' ') {
565 if (mc_ != 0) {
566 std::copy(ma_, ma_ + s - 1, &m_ab[0]);
567 std::copy(mb_, mb_ + s - 1, &m_ab[s - 1]);
568 PARDISO_LU_seq_sparse_direct_solver<double>::resolve(
569 s - 1, 2, &m_ab[0], s - 1, &m_ab[0], s - 1);
570 div = 1 / (*mc_ - blas::dot(s - 1, ma_, 1, &m_ab[s - 1], 1));
571 }
572
573 for (size_t i = 0; i != nrhs; ++i) {
574 const double x2 = x[ldx * i + s - 1] =
575 div * (b[ldb * i + s - 1] - blas::dot(s - 1, b + ldb * i, 1, &m_ab[s - 1], 1));
576 PARDISO_LU_seq_sparse_direct_solver<double>::resolve(
577 s - 1, 1, b + ldb * i, s - 1, x + ldx * i, s - 1);
578 blas::axpy(s - 1, -x2, &m_ab[0], 1, x + ldx * i, 1);
579 }
580 }
581
582protected:
583 std::vector<double> m_ab;
584 double div;
585};
586
587}} // namespace b2000::b2linalg
588
589#endif /* B2PARDISO_SOLVER_H_ */
THROW
#define THROW
Definition b2exception.H:198
b2000::logging::get_logger
Logger & get_logger(const std::string &logger_name="")
Definition b2logging.H:829
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::UnimplementedError
GenericException< UnimplementedError_name > UnimplementedError
Definition b2exception.H:314