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