Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
XtrMumps.c
Go to the documentation of this file.
1
14#include <time.h>
15
16#include "fasp.h"
17#include "fasp_functs.h"
18
19#if WITH_MUMPS
20#include "dmumps_c.h"
21#endif
22
23#define ICNTL(I) icntl[(I)-1]
25/*---------------------------------*/
26/*-- Public Functions --*/
27/*---------------------------------*/
28
45int fasp_solver_mumps(dCSRmat* ptrA, dvector* b, dvector* u, const SHORT prtlvl)
46{
47
48#if WITH_MUMPS
49
50 DMUMPS_STRUC_C id;
51
52 int status = FASP_SUCCESS;
53 const int n = ptrA->row;
54 const int nz = ptrA->nnz;
55 int* IA = ptrA->IA;
56 int* JA = ptrA->JA;
57 double* AA = ptrA->val;
58 double* f = b->val;
59 double* x = u->val;
60
61 int* irn;
62 int* jcn;
63 double* a;
64 double* rhs;
65 int i, j;
66 int begin_row, end_row;
67
68#if DEBUG_MODE
69 printf("### DEBUG: fasp_solver_mumps ... [Start]\n");
70 printf("### DEBUG: nr=%d, nnz=%d\n", n, nz);
71#endif
72
73 // First check the matrix format
74 if (IA[0] != 0 && IA[0] != 1) {
75 printf("### ERROR: Matrix format is wrong -- IA[0] = %d\n", IA[0]);
76 return ERROR_SOLVER_EXIT;
77 }
78
79 REAL start_time, end_time;
80 fasp_gettime(&start_time);
81
82 /* Define A and rhs */
83 irn = (int*)malloc(sizeof(int) * nz);
84 jcn = (int*)malloc(sizeof(int) * nz);
85 a = (double*)malloc(sizeof(double) * nz);
86 rhs = (double*)malloc(sizeof(double) * n);
87
88 if (IA[0] == 0) { // C-convention
89 for (i = 0; i < n; i++) {
90 begin_row = IA[i];
91 end_row = IA[i + 1];
92 for (j = begin_row; j < end_row; j++) {
93 irn[j] = i + 1;
94 jcn[j] = JA[j] + 1;
95 a[j] = AA[j];
96 }
97 }
98 } else { // F-convention
99 for (i = 0; i < n; i++) {
100 begin_row = IA[i] - 1;
101 end_row = IA[i + 1] - 1;
102 for (j = begin_row; j < end_row; j++) {
103 irn[j] = i + 1;
104 jcn[j] = JA[j];
105 a[j] = AA[j];
106 }
107 }
108 }
109
110 /* Initialize a MUMPS instance. */
111 id.job = -1;
112 id.par = 1; // host involved in factorization/solve
113 id.sym = 0; // 0: general, 1: spd, 2: sym
114 id.comm_fortran = 0;
115 dmumps_c(&id);
116
117 /* Define the problem on the host */
118 id.n = n;
119 id.nz = nz;
120 id.irn = irn;
121 id.jcn = jcn;
122 id.a = a;
123 id.rhs = rhs;
124
125 if (prtlvl > PRINT_MORE) { // enable debug output
126 id.ICNTL(1) = 6; // err output stream
127 id.ICNTL(2) = 6; // warn/info output stream
128 id.ICNTL(3) = 6; // global output stream
129 id.ICNTL(4) = 3; // 0:none, 1:err, 2:warn/stats, 3:diagnos, 4:parameters
130 if (prtlvl > PRINT_MOST) id.ICNTL(11) = 1; // 0:none, 1:all, 2:main
131 } else {
132 id.ICNTL(1) = -1;
133 id.ICNTL(2) = -1;
134 id.ICNTL(3) = -1;
135 id.ICNTL(4) = 0;
136 }
137
138 /* Call the MUMPS package */
139 for (i = 0; i < n; i++) rhs[i] = f[i];
140
141 id.job = 6; /* Combines phase 1, 2, and 3 */
142 dmumps_c(&id); /* Sometimes segmentation faults in MUMPS-5.0.0 */
143 status = id.info[0];
144 if (status < 0) {
145 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
146 printf("### ERROR: MUMPS solve failed!\n");
147 exit(ERROR_SOLVER_MISC);
148 }
149
150 for (i = 0; i < n; i++) x[i] = id.rhs[i];
151
152 id.job = -2;
153 dmumps_c(&id); /* Terminate instance */
154
155 free(irn);
156 free(jcn);
157 free(a);
158 free(rhs);
159
160 if (prtlvl > PRINT_MIN) {
161 fasp_gettime(&end_time);
162 fasp_cputime("MUMPS solver", end_time - start_time);
163 }
164
165#if DEBUG_MODE
166 printf("### DEBUG: fasp_solver_mumps ... [Finish]\n");
167#endif
168 return status;
169
170#else
171
172 printf("### ERROR: MUMPS is not available!\n");
173 return ERROR_SOLVER_EXIT;
174
175#endif
176}
177
196int fasp_solver_mumps_steps(dCSRmat* ptrA, dvector* b, dvector* u, Mumps_data* mumps)
197{
198#if WITH_MUMPS
199
200 DMUMPS_STRUC_C id;
201
202 int status = FASP_SUCCESS;
203 int job = mumps->job;
204 static int job_stat = 0;
205 int i, j;
206 int* irn;
207 int* jcn;
208 double* a;
209 double* rhs;
210
211 switch (job) {
212
213 case 1:
214 {
215#if DEBUG_MODE
216 printf("### DEBUG: %s, step %d, job_stat = %d... [Start]\n",
217 __FUNCTION__, job, job_stat);
218#endif
219 int begin_row, end_row;
220 const int n = ptrA->row;
221 const int nz = ptrA->nnz;
222 int* IA = ptrA->IA;
223 int* JA = ptrA->JA;
224 double* AA = ptrA->val;
225
226 irn = id.irn = (int*)malloc(sizeof(int) * nz);
227 jcn = id.jcn = (int*)malloc(sizeof(int) * nz);
228 a = id.a = (double*)malloc(sizeof(double) * nz);
229 rhs = id.rhs = (double*)malloc(sizeof(double) * n);
230
231 // First check the matrix format
232 if (IA[0] != 0 && IA[0] != 1) {
233 printf("### ERROR: Matrix format is wrong, IA[0] = %d!\n", IA[0]);
234 return ERROR_SOLVER_EXIT;
235 }
236
237 // Define A and rhs
238 if (IA[0] == 0) { // C-convention
239 for (i = 0; i < n; i++) {
240 begin_row = IA[i];
241 end_row = IA[i + 1];
242 for (j = begin_row; j < end_row; j++) {
243 irn[j] = i + 1;
244 jcn[j] = JA[j] + 1;
245 a[j] = AA[j];
246 }
247 }
248 } else { // F-convention
249 for (i = 0; i < n; i++) {
250 begin_row = IA[i] - 1;
251 end_row = IA[i + 1] - 1;
252 for (j = begin_row; j < end_row; j++) {
253 irn[j] = i + 1;
254 jcn[j] = JA[j];
255 a[j] = AA[j];
256 }
257 }
258 }
259
260 /* Initialize a MUMPS instance. */
261 id.job = -1;
262 id.par = 1;
263 id.sym = 0;
264 id.comm_fortran = 0;
265 dmumps_c(&id);
266
267 /* Define the problem on the host */
268 id.n = n;
269 id.nz = nz;
270 id.irn = irn;
271 id.jcn = jcn;
272 id.a = a;
273 id.rhs = rhs;
274
275 /* No outputs */
276 id.ICNTL(1) = -1;
277 id.ICNTL(2) = -1;
278 id.ICNTL(3) = -1;
279 id.ICNTL(4) = 0;
280
281 id.job = 4;
282 dmumps_c(&id);
283 status = id.info[0];
284 if (status < 0) {
285 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
286 printf("### ERROR: MUMPS factorization failed!\n");
287 exit(ERROR_SOLVER_MISC);
288 }
289 job_stat = 1;
290
291 mumps->id = id;
292
293#if DEBUG_MODE
294 printf("### DEBUG: %s, step %d, job_stat = %d... [Finish]\n",
295 __FUNCTION__, job, job_stat);
296#endif
297 break;
298 }
299
300 case 2:
301 {
302#if DEBUG_MODE
303 printf("### DEBUG: %s, step %d, job_stat = %d... [Start]\n",
304 __FUNCTION__, job, job_stat);
305#endif
306 id = mumps->id;
307
308 if (job_stat != 1)
309 printf("### ERROR: %s setup failed!\n", __FUNCTION__);
310
311 /* Call the MUMPS package. */
312 for (i = 0; i < id.n; i++) id.rhs[i] = b->val[i];
313
314 id.job = 3;
315 dmumps_c(&id);
316 status = id.info[0];
317 if (status < 0) {
318 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
319 printf("### ERROR: MUMPS solve failed!\n");
320 exit(ERROR_SOLVER_MISC);
321 }
322
323 for (i = 0; i < id.n; i++) u->val[i] = id.rhs[i];
324
325#if DEBUG_MODE
326 printf("### DEBUG: %s, step %d, job_stat = %d... [Finish]\n",
327 __FUNCTION__, job, job_stat);
328#endif
329 break;
330 }
331
332 case 3:
333 {
334#if DEBUG_MODE
335 printf("### DEBUG: %s, step %d, job_stat = %d... [Start]\n",
336 __FUNCTION__, job, job_stat);
337#endif
338 id = mumps->id;
339
340 if (job_stat != 1)
341 printf("### ERROR: %s setup failed!\n", __FUNCTION__);
342
343 free(id.irn);
344 free(id.jcn);
345 free(id.a);
346 free(id.rhs);
347 id.job = -2;
348 dmumps_c(&id); /* Terminate instance */
349
350#if DEBUG_MODE
351 printf("### DEBUG: %s, step %d, job_stat = %d... [Finish]\n",
352 __FUNCTION__, job, job_stat);
353#endif
354
355 break;
356 }
357
358 default:
359 printf("### ERROR: job = %d. Should be 1, 2, or 3!\n", job);
360 return ERROR_SOLVER_EXIT;
361 }
362
363 return status;
364
365#else
366
367 printf("### ERROR: MUMPS is not available!\n");
368 return ERROR_SOLVER_EXIT;
369
370#endif
371}
372
373#if WITH_MUMPS
387Mumps_data fasp_mumps_factorize(dCSRmat* ptrA, dvector* b, dvector* u,
388 const SHORT prtlvl)
389{
390 Mumps_data mumps;
391 DMUMPS_STRUC_C id;
392
393 int status = FASP_SUCCESS;
394 const int n = ptrA->col;
395 const int nz = ptrA->nnz;
396 int* IA = ptrA->IA;
397 int* JA = ptrA->JA;
398 double* AA = ptrA->val;
399 int i, j;
400
401 int* irn = id.irn = (int*)malloc(sizeof(int) * nz);
402 int* jcn = id.jcn = (int*)malloc(sizeof(int) * nz);
403 double* a = id.a = (double*)malloc(sizeof(double) * nz);
404 double* rhs = id.rhs = (double*)malloc(sizeof(double) * n);
405
406 int begin_row, end_row;
407
408#if DEBUG_MODE
409 printf("### DEBUG: %s ... [Start]\n", __FUNCTION__);
410 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", n, n, nz);
411#endif
412
413 clock_t start_time = clock();
414
415 if (IA[0] == 0) { // C-convention
416 for (i = 0; i < n; i++) {
417 begin_row = IA[i];
418 end_row = IA[i + 1];
419 for (j = begin_row; j < end_row; j++) {
420 irn[j] = i + 1;
421 jcn[j] = JA[j] + 1;
422 a[j] = AA[j];
423 }
424 }
425 } else { // F-convention
426 for (i = 0; i < n; i++) {
427 begin_row = IA[i] - 1;
428 end_row = IA[i + 1] - 1;
429 for (j = begin_row; j < end_row; j++) {
430 irn[j] = i + 1;
431 jcn[j] = JA[j];
432 a[j] = AA[j];
433 }
434 }
435 }
436
437 /* Initialize a MUMPS instance. */
438 id.job = -1;
439 id.par = 1;
440 id.sym = 0;
441 id.comm_fortran = 0;
442 dmumps_c(&id);
443
444 /* Define the problem on the host */
445 id.n = n;
446 id.nz = nz;
447 id.irn = irn;
448 id.jcn = jcn;
449 id.a = a;
450 id.rhs = rhs;
451
452 if (prtlvl > PRINT_MORE) { // enable debug output
453 id.ICNTL(1) = 6; // err output stream
454 id.ICNTL(2) = 6; // warn/info output stream
455 id.ICNTL(3) = 6; // global output stream
456 id.ICNTL(4) = 3; // 0:none, 1:err, 2:warn/stats, 3:diagnos, 4:parameters
457 if (prtlvl > PRINT_MOST) id.ICNTL(11) = 1; // 0:none, 1:all, 2:main
458 } else {
459 id.ICNTL(1) = -1;
460 id.ICNTL(2) = -1;
461 id.ICNTL(3) = -1;
462 id.ICNTL(4) = 0;
463 }
464
465 id.job = 4;
466 dmumps_c(&id);
467 status = id.info[0];
468 if (status < 0) {
469 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
470 printf("### ERROR: MUMPS factorization failed!\n");
471 exit(ERROR_SOLVER_MISC);
472 }
473
474 if (prtlvl > PRINT_MIN) {
475 clock_t end_time = clock();
476 double fac_time = (double)(end_time - start_time) / (double)(CLOCKS_PER_SEC);
477 printf("MUMPS factorize costs %f seconds.\n", fac_time);
478 }
479
480#if DEBUG_MODE
481 printf("### DEBUG: %s ... [Finish]\n", __FUNCTION__);
482#endif
483
484 mumps.id = id;
485
486 return mumps;
487}
488#endif
489
490#if WITH_MUMPS
505void fasp_mumps_solve(dCSRmat* ptrA, dvector* b, dvector* u, Mumps_data mumps,
506 const SHORT prtlvl)
507{
508 int i;
509 DMUMPS_STRUC_C id = mumps.id;
510 int status = FASP_SUCCESS;
511 double* rhs = id.rhs;
512
513#if DEBUG_MODE
514 printf("### DEBUG: %s ... [Start]\n", __FUNCTION__);
515 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nz);
516#endif
517
518 clock_t start_time = clock();
519
520 double* f = b->val;
521 double* x = u->val;
522
523 /* Call the MUMPS package. */
524 for (i = 0; i < id.n; i++) rhs[i] = f[i];
525
526 if (prtlvl > PRINT_MORE) { // enable debug output
527 id.ICNTL(1) = 6; // err output stream
528 id.ICNTL(2) = 6; // warn/info output stream
529 id.ICNTL(3) = 6; // global output stream
530 id.ICNTL(4) = 3; // 0:none, 1:err, 2:warn/stats, 3:diagnos, 4:parameters
531 if (prtlvl > PRINT_MOST) id.ICNTL(11) = 1; // 0:none, 1:all, 2:main
532 } else {
533 id.ICNTL(1) = -1;
534 id.ICNTL(2) = -1;
535 id.ICNTL(3) = -1;
536 id.ICNTL(4) = 0;
537 }
538
539 id.job = 3;
540 dmumps_c(&id);
541 status = id.info[0];
542 if (status < 0) {
543 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
544 printf("### ERROR: MUMPS solve failed!\n");
545 exit(ERROR_SOLVER_MISC);
546 }
547
548 for (i = 0; i < id.n; i++) x[i] = id.rhs[i];
549
550 if (prtlvl > PRINT_NONE) {
551 clock_t end_time = clock();
552 double solve_time = (double)(end_time - start_time) / (double)(CLOCKS_PER_SEC);
553 printf("MUMPS costs %f seconds.\n", solve_time);
554 }
555
556#if DEBUG_MODE
557 printf("### DEBUG: %s ... [Finish]\n", __FUNCTION__);
558#endif
559}
560#endif
561
562#if WITH_MUMPS
573void fasp_mumps_free(Mumps_data* mumps)
574{
575 DMUMPS_STRUC_C id = mumps->id;
576
577 free(id.irn);
578 free(id.jcn);
579 free(id.a);
580 free(id.rhs);
581}
582#endif
583
584/*---------------------------------*/
585/*-- End of File --*/
586/*---------------------------------*/
fasp_cputime
void fasp_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: AuxMessage.c:179
fasp_gettime
void fasp_gettime(REAL *time)
Get system time.
Definition: AuxTiming.c:37
fasp_solver_mumps
int fasp_solver_mumps(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Ax=b by MUMPS directly.
Definition: XtrMumps.c:45
fasp_solver_mumps_steps
int fasp_solver_mumps_steps(dCSRmat *ptrA, dvector *b, dvector *u, Mumps_data *mumps)
Solve Ax=b by MUMPS in three steps.
Definition: XtrMumps.c:196
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
PRINT_MOST
#define PRINT_MOST
Definition: fasp_const.h:77
FASP_SUCCESS
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
ERROR_SOLVER_MISC
#define ERROR_SOLVER_MISC
Definition: fasp_const.h:47
PRINT_MORE
#define PRINT_MORE
Definition: fasp_const.h:76
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
ERROR_SOLVER_EXIT
#define ERROR_SOLVER_EXIT
Definition: fasp_const.h:49
PRINT_MIN
#define PRINT_MIN
Definition: fasp_const.h:74
Mumps_data
Data for MUMPS interface.
Definition: fasp.h:607
Mumps_data::job
INT job
work for MUMPS
Definition: fasp.h:615
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:163
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:166
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360