Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
interface_mumps.c
Go to the documentation of this file.
1 
8 #include <time.h>
9 
10 #include "fasp.h"
11 #include "fasp_functs.h"
12 
13 #if WITH_MUMPS
14 #include "dmumps_c.h"
15 #endif
16 
17 #define ICNTL(I) icntl[(I)-1]
19 /*---------------------------------*/
20 /*-- Public Functions --*/
21 /*---------------------------------*/
22 
40  dvector *b,
41  dvector *u,
42  const SHORT prtlvl)
43 {
44 
45 #if WITH_MUMPS
46 
47  DMUMPS_STRUC_C id;
48 
49  const int n = ptrA->row;
50  const int nz = ptrA->nnz;
51  int *IA = ptrA->IA;
52  int *JA = ptrA->JA;
53  double *AA = ptrA->val;
54  double *b1 = b->val;
55  double *x = u->val;
56 
57  int *irn;
58  int *jcn;
59  double *a;
60  double *rhs;
61  int i,j;
62  int begin_row, end_row;
63 
64 #if DEBUG_MODE
65  printf("### DEBUG: fasp_solver_mumps ... [Start]\n");
66  printf("### DEBUG: nr=%d, nnz=%d\n", n, nz);
67 #endif
68 
69  // First check the matrix format
70  if ( IA[0] != 0 && IA[0] != 1 ) {
71  printf("### ERROR: Matrix format is wrong -- IA[0] = %d\n", IA[0]);
72  return ERROR_SOLVER_EXIT;
73  }
74 
75  clock_t start_time = clock();
76 
77  /* Define A and rhs */
78  irn = (int *)malloc( sizeof(int)*nz );
79  jcn = (int *)malloc( sizeof(int)*nz );
80  a = (double *)malloc( sizeof(double)*nz );
81  rhs = (double *)malloc( sizeof(double)*n );
82 
83  if ( IA[0] == 0 ) { // C-convention
84  for (i=0; i<n; i++) {
85  begin_row = IA[i]; end_row = IA[i+1];
86  for (j=begin_row; j< end_row; j++) {
87  irn[j] = i + 1;
88  jcn[j] = JA[j]+1;
89  a[j] = AA[j];
90  }
91  }
92  }
93  else { // For-convention
94  for (i=0; i<n; i++) {
95  begin_row = IA[i]-1; end_row = IA[i+1]-1;
96  for (j=begin_row; j< end_row; j++) {
97  irn[j] = i + 1;
98  jcn[j] = JA[j];
99  a[j] = AA[j];
100  }
101  }
102  }
103 
104  /* Initialize a MUMPS instance. */
105  id.job=-1; id.par=1; id.sym=0; id.comm_fortran=0;
106  dmumps_c(&id);
107 
108  /* Define the problem on the host */
109  id.n = n; id.nz =nz; id.irn=irn; id.jcn=jcn;
110  id.a = a; id.rhs = rhs;
111 
112  /* No outputs */
113  id.ICNTL(1) = -1;
114  id.ICNTL(2) = -1;
115  id.ICNTL(3) = -1;
116  id.ICNTL(4) = 0;
117 
118  /* Call the MUMPS package */
119  for (i=0; i<n; i++) rhs[i] = b1[i];
120 
121  id.job = 6; dmumps_c(&id);
122 
123  for (i=0; i<n; i++) x[i] = id.rhs[i];
124 
125  id.job = -2;
126  dmumps_c(&id); /* Terminate instance */
127 
128  free(irn);
129  free(jcn);
130  free(a);
131  free(rhs);
132 
133  if ( prtlvl > PRINT_MIN ) {
134  clock_t end_time = clock();
135  double solve_time = (double)(end_time - start_time)/(double)(CLOCKS_PER_SEC);
136  printf("MUMPS costs %f seconds.\n", solve_time);
137  }
138 
139 #if DEBUG_MODE
140  printf("### DEBUG: fasp_solver_mumps ... [Finish]\n");
141 #endif
142  return FASP_SUCCESS;
143 #else
144 
145  printf("### ERROR: MUMPS is not available!\n");
146  return ERROR_SOLVER_EXIT;
147 
148 #endif
149 
150 }
151 
170  dvector *b,
171  dvector *u,
172  Mumps_data *mumps)
173 {
174 #if WITH_MUMPS
175 
176  DMUMPS_STRUC_C id;
177 
178  int job = mumps->job;
179 
180  static int job_stat = 0;
181  int i,j;
182 
183  int *irn;
184  int *jcn;
185  double *a;
186  double *rhs;
187 
188 #if DEBUG_MODE
189  printf("### DEBUG: %s job_stat = %d\n", __FUNCTION__, job_stat);
190 #endif
191 
192  switch ( job ) {
193 
194  case 1:
195  {
196 #if DEBUG_MODE
197  printf("### DEBUG: %s Step 1 ... [Start]\n", __FUNCTION__);
198 #endif
199  int begin_row, end_row;
200  const int n = ptrA->row;
201  const int nz = ptrA->nnz;
202  int *IA = ptrA->IA;
203  int *JA = ptrA->JA;
204  double *AA = ptrA->val;
205 
206  irn = id.irn = (int *)malloc( sizeof(int)*nz );
207  jcn = id.jcn = (int *)malloc( sizeof(int)*nz );
208  a = id.a = (double *)malloc( sizeof(double)*nz );
209  rhs = id.rhs = (double *)malloc( sizeof(double)*n );
210 
211  // First check the matrix format
212  if ( IA[0] != 0 && IA[0] != 1 ) {
213  printf("### ERROR: Matrix format is wrong -- IA[0] = %d\n", IA[0]);
214  return ERROR_SOLVER_EXIT;
215  }
216 
217  // Define A and rhs
218  if ( IA[0] == 0 ) { // C-convention
219  for (i=0; i<n; i++) {
220  begin_row = IA[i]; end_row = IA[i+1];
221  for (j=begin_row; j< end_row; j++) {
222  irn[j] = i + 1;
223  jcn[j] = JA[j]+1;
224  a[j] = AA[j];
225  }
226  }
227  }
228  else { // For-convention
229  for (i=0; i<n; i++) {
230  begin_row = IA[i]-1; end_row = IA[i+1]-1;
231  for (j=begin_row; j< end_row; j++) {
232  irn[j] = i + 1;
233  jcn[j] = JA[j];
234  a[j] = AA[j];
235  }
236  }
237  }
238 
239  /* Initialize a MUMPS instance. */
240  id.job = -1; id.par = 1; id.sym = 0; id.comm_fortran = 0;
241  dmumps_c(&id);
242  /* Define the problem on the host */
243  id.n = n; id.nz = nz; id.irn = irn; id.jcn = jcn;
244  id.a = a; id.rhs = rhs;
245 
246  /* No outputs */
247  id.ICNTL(1) = -1;
248  id.ICNTL(2) = -1;
249  id.ICNTL(3) = -1;
250  id.ICNTL(4) = 0;
251 
252  id.job = 4; dmumps_c(&id);
253  job_stat = 1;
254 
255  mumps->id = id;
256 
257 #if DEBUG_MODE
258  printf("### DEBUG: %s, Step 1 ... [Finish]\n", __FUNCTION__);
259 #endif
260  break;
261  }
262 
263  case 2:
264  {
265 #if DEBUG_MODE
266  printf("### DEBUG: %s, Step 2 ... [Start]\n", __FUNCTION__);
267 #endif
268  id = mumps->id;
269 
270  if ( job_stat != 1 )
271  printf("### ERROR: %s setup not finished!\n", __FUNCTION__);
272 
273  /* Call the MUMPS package. */
274  for(i=0; i<id.n; i++) id.rhs[i] = b->val[i];
275 
276  id.job=3; dmumps_c(&id);
277 
278  for(i=0; i<id.n; i++) u->val[i] = id.rhs[i];
279 
280 #if DEBUG_MODE
281  printf("### DEBUG: %s, Step 2 ... [Finish]\n", __FUNCTION__);
282 #endif
283  break;
284  }
285 
286  case 3:
287  {
288  id = mumps->id;
289 
290  if ( job_stat !=1 )
291  printf("### ERROR: %s setup not finished!\n", __FUNCTION__);
292 
293  free(id.irn);
294  free(id.jcn);
295  free(id.a);
296  free(id.rhs);
297  id.job = -2;
298  dmumps_c(&id); /* Terminate instance */
299 
300  break;
301  }
302 
303  default:
304  printf("### ERROR: Parameter job must be 1, 2, or 3!\n");
305  return ERROR_SOLVER_EXIT;
306 
307  }
308 
309  return FASP_SUCCESS;
310 
311 #else
312 
313  printf("### ERROR: MUMPS is not available!\n");
314  return ERROR_SOLVER_EXIT;
315 
316 #endif
317 }
318 
319 #if WITH_MUMPS
320 
333 Mumps_data fasp_mumps_factorize (dCSRmat *ptrA,
334  dvector *b,
335  dvector *u,
336  const SHORT prtlvl)
337 {
338  Mumps_data mumps;
339  DMUMPS_STRUC_C id;
340 
341  int i,j;
342  const int m = ptrA->row;
343  const int n = ptrA->col;
344  const int nz = ptrA->nnz;
345  int *IA = ptrA->IA;
346  int *JA = ptrA->JA;
347  double *AA = ptrA->val;
348 
349  int *irn = id.irn = (int *)malloc( sizeof(int)*nz );
350  int *jcn = id.jcn = (int *)malloc( sizeof(int)*nz );
351  double *a = id.a = (double *)malloc( sizeof(double)*nz );
352  double *rhs = id.rhs = (double *)malloc( sizeof(double)*n );
353 
354  int begin_row, end_row;
355 
356 #if DEBUG_MODE
357  printf("### DEBUG: %s ... [Start]\n", __FUNCTION__);
358  printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nz);
359 #endif
360 
361  clock_t start_time = clock();
362 
363  if ( IA[0] == 0 ) { // C-convention
364  for (i=0; i<n; i++) {
365  begin_row = IA[i]; end_row = IA[i+1];
366  for (j=begin_row; j< end_row; j++) {
367  irn[j] = i + 1;
368  jcn[j] = JA[j]+1;
369  a[j] = AA[j];
370  }
371  }
372  }
373  else { // For-convention
374  for (i=0; i<n; i++) {
375  begin_row = IA[i]-1; end_row = IA[i+1]-1;
376  for (j=begin_row; j< end_row; j++) {
377  irn[j] = i + 1;
378  jcn[j] = JA[j];
379  a[j] = AA[j];
380  }
381  }
382  }
383 
384  /* Initialize a MUMPS instance. */
385  id.job = -1; id.par = 1; id.sym = 0; id.comm_fortran = 0;
386  dmumps_c(&id);
387  /* Define the problem on the host */
388  id.n = n; id.nz = nz; id.irn = irn; id.jcn = jcn;
389  id.a = a; id.rhs = rhs;
390 
391  /* No outputs */
392  id.ICNTL(1) = -1;
393  id.ICNTL(2) = -1;
394  id.ICNTL(3) = -1;
395  id.ICNTL(4) = 0;
396 
397  id.job = 4; dmumps_c(&id);
398 
399  if ( prtlvl > PRINT_MIN ) {
400  clock_t end_time = clock();
401  double fac_time = (double)(end_time - start_time)/(double)(CLOCKS_PER_SEC);
402  printf("UMFPACK factorize costs %f seconds.\n", fac_time);
403  }
404 
405 #if DEBUG_MODE
406  printf("### DEBUG: %s ... [Finish]\n", __FUNCTION__);
407 #endif
408 
409  mumps.id = id;
410 
411  return mumps;
412 }
413 #endif
414 
415 #if WITH_MUMPS
416 
430 void fasp_mumps_solve (dCSRmat *ptrA,
431  dvector *b,
432  dvector *u,
433  Mumps_data mumps,
434  const SHORT prtlvl)
435 {
436  int i,j;
437 
438  DMUMPS_STRUC_C id = mumps.id;
439 
440  const int m = ptrA->row;
441  const int n = ptrA->row;
442  const int nz = ptrA->nnz;
443  int *IA = ptrA->IA;
444  int *JA = ptrA->JA;
445  double *AA = ptrA->val;
446 
447  int *irn = id.irn;
448  int *jcn = id.jcn;
449  double *a = id.a;
450  double *rhs = id.rhs;
451 
452 #if DEBUG_MODE
453  printf("### DEBUG: %s ... [Start]\n", __FUNCTION__);
454  printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nz);
455 #endif
456 
457  clock_t start_time = clock();
458 
459  double *b1 = b->val;
460  double *x = u->val;
461 
462  /* Call the MUMPS package. */
463  for(i=0; i<id.n; i++) rhs[i] = b1[i];
464 
465  id.job=3; dmumps_c(&id);
466 
467  for(i=0; i<id.n; i++) x[i] = id.rhs[i];
468 
469  if ( prtlvl > PRINT_NONE ) {
470  clock_t end_time = clock();
471  double solve_time = (double)(end_time - start_time)/(double)(CLOCKS_PER_SEC);
472  printf("UMFPACK costs %f seconds.\n", solve_time);
473  }
474 
475 #if DEBUG_MODE
476  printf("### DEBUG: %s ... [Finish]\n", __FUNCTION__);
477 #endif
478 }
479 #endif
480 
481 #if WITH_MUMPS
482 
492 void fasp_mumps_free (Mumps_data *mumps)
493 {
494  DMUMPS_STRUC_C id = mumps->id;
495 
496  free(id.irn);
497  free(id.jcn);
498  free(id.a);
499  free(id.rhs);
500 }
501 #endif
502 
503 /*---------------------------------*/
504 /*-- End of File --*/
505 /*---------------------------------*/
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
int fasp_solver_mumps_steps(dCSRmat *ptrA, dvector *b, dvector *u, Mumps_data *mumps)
Solve Ax=b by MUMPS in three steps.
REAL * val
nonzero entries of A
Definition: fasp.h:166
REAL * val
actual vector entries
Definition: fasp.h:348
Vector with n entries of REAL type.
Definition: fasp.h:342
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
Parameters for MUMPS interface.
Definition: fasp.h:459
INT nnz
number of nonzero entries
Definition: fasp.h:157
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
int fasp_solver_mumps(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Ax=b by MUMPS directly.
#define ERROR_SOLVER_EXIT
Definition: fasp_const.h:55
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
INT job
work for MUMPS
Definition: fasp.h:467
#define PRINT_MIN
Definition: fasp_const.h:80