xref: /petsc/src/mat/impls/aij/mpi/mumps/mumps.c (revision d32f9abdbc052d6e1fd06679b17a55415c3aae30) 
1 #define PETSCMAT_DLL
2 
3 /*
4     Provides an interface to the MUMPS sparse solver
5 */
6 #include "src/mat/impls/aij/seq/aij.h"
7 #include "src/mat/impls/aij/mpi/mpiaij.h"
8 #include "src/mat/impls/sbaij/seq/sbaij.h"
9 #include "src/mat/impls/sbaij/mpi/mpisbaij.h"
10 
11 EXTERN_C_BEGIN
12 #if defined(PETSC_USE_COMPLEX)
13 #include "zmumps_c.h"
14 #else
15 #include "dmumps_c.h"
16 #endif
17 EXTERN_C_END
18 #define JOB_INIT -1
19 #define JOB_END -2
20 /* macros s.t. indices match MUMPS documentation */
21 #define ICNTL(I) icntl[(I)-1]
22 #define CNTL(I) cntl[(I)-1]
23 #define INFOG(I) infog[(I)-1]
24 #define INFO(I) info[(I)-1]
25 #define RINFOG(I) rinfog[(I)-1]
26 #define RINFO(I) rinfo[(I)-1]
27 
28 typedef struct {
29 #if defined(PETSC_USE_COMPLEX)
30   ZMUMPS_STRUC_C id;
31 #else
32   DMUMPS_STRUC_C id;
33 #endif
34   MatStructure   matstruc;
35   PetscMPIInt    myid,size;
36   PetscInt       *irn,*jcn,sym,nSolve;
37   PetscScalar    *val;
38   MPI_Comm       comm_mumps;
39   VecScatter     scat_rhs, scat_sol;
40   PetscTruth     isAIJ,CleanUpMUMPS;
41   Vec            b_seq,x_seq;
42   PetscErrorCode (*MatDestroy)(Mat);
43 } Mat_MUMPS;
44 
45 EXTERN PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*);
46 
47 /* convert Petsc mpiaij matrix to triples: row[nz], col[nz], val[nz] */
48 /*
49   input:
50     A       - matrix in mpiaij or mpisbaij (bs=1) format
51     shift   - 0: C style output triple; 1: Fortran style output triple.
52     valOnly - FALSE: spaces are allocated and values are set for the triple
53               TRUE:  only the values in v array are updated
54   output:
55     nnz     - dim of r, c, and v (number of local nonzero entries of A)
56     r, c, v - row and col index, matrix values (matrix triples)
57  */
58 PetscErrorCode MatConvertToTriples(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v)
59 {
60   PetscInt       *ai, *aj, *bi, *bj, rstart,nz, *garray;
61   PetscErrorCode ierr;
62   PetscInt       i,j,jj,jB,irow,m=A->rmap.n,*ajj,*bjj,countA,countB,colA_start,jcol;
63   PetscInt       *row,*col;
64   PetscScalar    *av, *bv,*val;
65   PetscTruth     isAIJ;
66 
67   PetscFunctionBegin;
68   ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isAIJ);CHKERRQ(ierr);
69   if (isAIJ){
70     Mat_MPIAIJ    *mat =  (Mat_MPIAIJ*)A->data;
71     Mat_SeqAIJ    *aa=(Mat_SeqAIJ*)(mat->A)->data;
72     Mat_SeqAIJ    *bb=(Mat_SeqAIJ*)(mat->B)->data;
73     nz = aa->nz + bb->nz;
74     ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap.rstart;
75     garray = mat->garray;
76     av=aa->a; bv=bb->a;
77 
78   } else {
79     Mat_MPISBAIJ  *mat =  (Mat_MPISBAIJ*)A->data;
80     Mat_SeqSBAIJ  *aa=(Mat_SeqSBAIJ*)(mat->A)->data;
81     Mat_SeqBAIJ    *bb=(Mat_SeqBAIJ*)(mat->B)->data;
82     if (A->rmap.bs > 1) SETERRQ1(PETSC_ERR_SUP," bs=%d is not supported yet\n", A->rmap.bs);
83     nz = aa->nz + bb->nz;
84     ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap.rstart;
85     garray = mat->garray;
86     av=aa->a; bv=bb->a;
87   }
88 
89   if (!valOnly){
90     ierr = PetscMalloc(nz*sizeof(PetscInt) ,&row);CHKERRQ(ierr);
91     ierr = PetscMalloc(nz*sizeof(PetscInt),&col);CHKERRQ(ierr);
92     ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr);
93     *r = row; *c = col; *v = val;
94   } else {
95     row = *r; col = *c; val = *v;
96   }
97   *nnz = nz;
98 
99   jj = 0; irow = rstart;
100   for ( i=0; i<m; i++ ) {
101     ajj = aj + ai[i];                 /* ptr to the beginning of this row */
102     countA = ai[i+1] - ai[i];
103     countB = bi[i+1] - bi[i];
104     bjj = bj + bi[i];
105 
106     /* get jB, the starting local col index for the 2nd B-part */
107     colA_start = rstart + ajj[0]; /* the smallest col index for A */
108     j=-1;
109     do {
110       j++;
111       if (j == countB) break;
112       jcol = garray[bjj[j]];
113     } while (jcol < colA_start);
114     jB = j;
115 
116     /* B-part, smaller col index */
117     colA_start = rstart + ajj[0]; /* the smallest col index for A */
118     for (j=0; j<jB; j++){
119       jcol = garray[bjj[j]];
120       if (!valOnly){
121         row[jj] = irow + shift; col[jj] = jcol + shift;
122 
123       }
124       val[jj++] = *bv++;
125     }
126     /* A-part */
127     for (j=0; j<countA; j++){
128       if (!valOnly){
129         row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift;
130       }
131       val[jj++] = *av++;
132     }
133     /* B-part, larger col index */
134     for (j=jB; j<countB; j++){
135       if (!valOnly){
136         row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift;
137       }
138       val[jj++] = *bv++;
139     }
140     irow++;
141   }
142 
143   PetscFunctionReturn(0);
144 }
145 
146 #undef __FUNCT__
147 #define __FUNCT__ "MatDestroy_MUMPS"
148 PetscErrorCode MatDestroy_MUMPS(Mat A)
149 {
150   Mat_MUMPS      *lu=(Mat_MUMPS*)A->spptr;
151   PetscErrorCode ierr;
152   PetscMPIInt    size=lu->size;
153 
154   PetscFunctionBegin;
155   if (lu->CleanUpMUMPS) {
156     /* Terminate instance, deallocate memories */
157     if (size > 1){
158       ierr = PetscFree(lu->id.sol_loc);CHKERRQ(ierr);
159       ierr = VecScatterDestroy(lu->scat_rhs);CHKERRQ(ierr);
160       ierr = VecDestroy(lu->b_seq);CHKERRQ(ierr);
161       if (lu->nSolve && lu->scat_sol){ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);}
162       if (lu->nSolve && lu->x_seq){ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);}
163       ierr = PetscFree(lu->val);CHKERRQ(ierr);
164     }
165     lu->id.job=JOB_END;
166 #if defined(PETSC_USE_COMPLEX)
167     zmumps_c(&lu->id);
168 #else
169     dmumps_c(&lu->id);
170 #endif
171     ierr = PetscFree(lu->irn);CHKERRQ(ierr);
172     ierr = PetscFree(lu->jcn);CHKERRQ(ierr);
173     ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr);
174   }
175   ierr = (lu->MatDestroy)(A);CHKERRQ(ierr);
176   PetscFunctionReturn(0);
177 }
178 
179 #undef __FUNCT__
180 #define __FUNCT__ "MatSolve_MUMPS"
181 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x)
182 {
183   Mat_MUMPS      *lu=(Mat_MUMPS*)A->spptr;
184   PetscScalar    *array;
185   Vec            x_seq;
186   IS             is_iden,is_petsc;
187   PetscErrorCode ierr;
188   PetscInt       i;
189 
190   PetscFunctionBegin;
191   lu->id.nrhs = 1;
192   x_seq = lu->b_seq;
193   if (lu->size > 1){
194     /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */
195     ierr = VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
196     ierr = VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
197     if (!lu->myid) {ierr = VecGetArray(x_seq,&array);CHKERRQ(ierr);}
198   } else {  /* size == 1 */
199     ierr = VecCopy(b,x);CHKERRQ(ierr);
200     ierr = VecGetArray(x,&array);CHKERRQ(ierr);
201   }
202   if (!lu->myid) { /* define rhs on the host */
203 #if defined(PETSC_USE_COMPLEX)
204     lu->id.rhs = (mumps_double_complex*)array;
205 #else
206     lu->id.rhs = array;
207 #endif
208   }
209   if (lu->size == 1){
210     ierr = VecRestoreArray(x,&array);CHKERRQ(ierr);
211   } else if (!lu->myid){
212     ierr = VecRestoreArray(x_seq,&array);CHKERRQ(ierr);
213   }
214 
215   if (lu->size > 1){
216     /* distributed solution */
217     lu->id.ICNTL(21) = 1;
218     if (!lu->nSolve){
219       /* Create x_seq=sol_loc for repeated use */
220       PetscInt    lsol_loc;
221       PetscScalar *sol_loc;
222       lsol_loc = lu->id.INFO(23); /* length of sol_loc */
223       ierr = PetscMalloc((1+lsol_loc)*(sizeof(PetscScalar)+sizeof(PetscInt)),&sol_loc);CHKERRQ(ierr);
224       lu->id.isol_loc = (PetscInt *)(sol_loc + lsol_loc);
225       lu->id.lsol_loc = lsol_loc;
226 #if defined(PETSC_USE_COMPLEX)
227       lu->id.sol_loc  = (ZMUMPS_DOUBLE *)sol_loc;
228 #else
229       lu->id.sol_loc  = (DMUMPS_DOUBLE *)sol_loc;
230 #endif
231       ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,lsol_loc,sol_loc,&lu->x_seq);CHKERRQ(ierr);
232     }
233   }
234 
235   /* solve phase */
236   /*-------------*/
237   lu->id.job = 3;
238 #if defined(PETSC_USE_COMPLEX)
239   zmumps_c(&lu->id);
240 #else
241   dmumps_c(&lu->id);
242 #endif
243   if (lu->id.INFOG(1) < 0) {
244     SETERRQ1(PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",lu->id.INFOG(1));
245   }
246 
247   if (lu->size > 1) { /* convert mumps distributed solution to petsc mpi x */
248     if (!lu->nSolve){ /* create scatter scat_sol */
249       ierr = ISCreateStride(PETSC_COMM_SELF,lu->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */
250       for (i=0; i<lu->id.lsol_loc; i++){
251         lu->id.isol_loc[i] -= 1; /* change Fortran style to C style */
252       }
253       ierr = ISCreateGeneral(PETSC_COMM_SELF,lu->id.lsol_loc,lu->id.isol_loc,&is_petsc);CHKERRQ(ierr);  /* to */
254       ierr = VecScatterCreate(lu->x_seq,is_iden,x,is_petsc,&lu->scat_sol);CHKERRQ(ierr);
255       ierr = ISDestroy(is_iden);CHKERRQ(ierr);
256       ierr = ISDestroy(is_petsc);CHKERRQ(ierr);
257     }
258     ierr = VecScatterBegin(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
259     ierr = VecScatterEnd(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
260   }
261   lu->nSolve++;
262   PetscFunctionReturn(0);
263 }
264 
265 #if !defined(PETSC_USE_COMPLEX)
266 /*
267   input:
268    F:        numeric factor
269   output:
270    nneg:     total number of negative pivots
271    nzero:    0
272    npos:     (global dimension of F) - nneg
273 */
274 
275 #undef __FUNCT__
276 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS"
277 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos)
278 {
279   Mat_MUMPS      *lu =(Mat_MUMPS*)F->spptr;
280   PetscErrorCode ierr;
281   PetscMPIInt    size;
282 
283   PetscFunctionBegin;
284   ierr = MPI_Comm_size(((PetscObject)F)->comm,&size);CHKERRQ(ierr);
285   /* MUMPS 4.3.1 calls ScaLAPACK when ICNTL(13)=0 (default), which does not offer the possibility to compute the inertia of a dense matrix. Set ICNTL(13)=1 to skip ScaLAPACK */
286   if (size > 1 && lu->id.ICNTL(13) != 1){
287     SETERRQ1(PETSC_ERR_ARG_WRONG,"ICNTL(13)=%d. -mat_mumps_icntl_13 must be set as 1 for correct global matrix inertia\n",lu->id.INFOG(13));
288   }
289   if (nneg){
290     if (!lu->myid){
291       *nneg = lu->id.INFOG(12);
292     }
293     ierr = MPI_Bcast(nneg,1,MPI_INT,0,lu->comm_mumps);CHKERRQ(ierr);
294   }
295   if (nzero) *nzero = 0;
296   if (npos)  *npos  = F->rmap.N - (*nneg);
297   PetscFunctionReturn(0);
298 }
299 #endif /* !defined(PETSC_USE_COMPLEX) */
300 
301 #undef __FUNCT__
302 #define __FUNCT__ "MatFactorNumeric_MUMPS"
303 PetscErrorCode MatFactorNumeric_MUMPS(Mat A,MatFactorInfo *info,Mat *F)
304 {
305   Mat_MUMPS      *lu =(Mat_MUMPS*)(*F)->spptr;
306   PetscErrorCode ierr;
307   PetscInt       rnz,nnz,nz=0,i,M=A->rmap.N,*ai,*aj,icntl;
308   PetscTruth     valOnly,flg;
309   Mat            F_diag;
310   IS             is_iden;
311   Vec            b;
312   PetscTruth     isSeqAIJ,isSeqSBAIJ;
313 
314   PetscFunctionBegin;
315   ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr);
316   ierr = PetscTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr);
317   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){
318     (*F)->ops->solve   = MatSolve_MUMPS;
319     (*F)->ops->destroy = MatDestroy_MUMPS;
320 
321     /* Initialize a MUMPS instance */
322     ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid);
323     ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr);
324     lu->id.job = JOB_INIT;
325     ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr);
326     lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps);
327 
328     /* Set mumps options */
329     ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr);
330     lu->id.par=1;  /* host participates factorizaton and solve */
331     lu->id.sym=lu->sym;
332     if (lu->sym == 2){
333       ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr);
334       if (flg && icntl == 1) lu->id.sym=icntl;  /* matrix is spd */
335     }
336 #if defined(PETSC_USE_COMPLEX)
337     zmumps_c(&lu->id);
338 #else
339     dmumps_c(&lu->id);
340 #endif
341 
342     if (isSeqAIJ || isSeqSBAIJ){
343       lu->id.ICNTL(18) = 0;   /* centralized assembled matrix input */
344     } else {
345       lu->id.ICNTL(18) = 3;   /* distributed assembled matrix input */
346     }
347 
348     icntl=-1;
349     lu->id.ICNTL(4) = 0;  /* level of printing; overwrite mumps default ICNTL(4)=2 */
350     ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr);
351     if ((flg && icntl > 0) || PetscLogPrintInfo) {
352       lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */
353     } else { /* no output */
354       lu->id.ICNTL(1) = 0;  /* error message, default= 6 */
355       lu->id.ICNTL(2) = -1; /* output stream for diagnostic printing, statistics, and warning. default=0 */
356       lu->id.ICNTL(3) = -1; /* output stream for global information, default=6 */
357     }
358     ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): matrix prescaling (0 to 7)","None",lu->id.ICNTL(6),&lu->id.ICNTL(6),PETSC_NULL);CHKERRQ(ierr);
359     icntl=-1;
360     ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr);
361     if (flg) {
362       if (icntl== 1){
363         SETERRQ(PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n");
364       } else {
365         lu->id.ICNTL(7) = icntl;
366       }
367     }
368     ierr = PetscOptionsInt("-mat_mumps_icntl_9","ICNTL(9): A or A^T x=b to be solved. 1: A; otherwise: A^T","None",lu->id.ICNTL(9),&lu->id.ICNTL(9),PETSC_NULL);CHKERRQ(ierr);
369     ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",lu->id.ICNTL(10),&lu->id.ICNTL(10),PETSC_NULL);CHKERRQ(ierr);
370     ierr = PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): error analysis, a positive value returns statistics (by -ksp_view)","None",lu->id.ICNTL(11),&lu->id.ICNTL(11),PETSC_NULL);CHKERRQ(ierr);
371     ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr);
372     ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control","None",lu->id.ICNTL(13),&lu->id.ICNTL(13),PETSC_NULL);CHKERRQ(ierr);
373     ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): percentage of estimated workspace increase","None",lu->id.ICNTL(14),&lu->id.ICNTL(14),PETSC_NULL);CHKERRQ(ierr);
374     ierr = PetscOptionsInt("-mat_mumps_icntl_15","ICNTL(15): efficiency control","None",lu->id.ICNTL(15),&lu->id.ICNTL(15),PETSC_NULL);CHKERRQ(ierr);
375 
376     ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr);
377     ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",lu->id.CNTL(2),&lu->id.CNTL(2),PETSC_NULL);CHKERRQ(ierr);
378     ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr);
379     ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",lu->id.CNTL(4),&lu->id.CNTL(4),PETSC_NULL);CHKERRQ(ierr);
380     PetscOptionsEnd();
381   }
382 
383   /* define matrix A */
384   switch (lu->id.ICNTL(18)){
385   case 0:  /* centralized assembled matrix input (size=1) */
386     if (!lu->myid) {
387       if (isSeqAIJ){
388         Mat_SeqAIJ   *aa = (Mat_SeqAIJ*)A->data;
389         nz               = aa->nz;
390         ai = aa->i; aj = aa->j; lu->val = aa->a;
391       } else if (isSeqSBAIJ) {
392         Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data;
393         nz                  =  aa->nz;
394         ai = aa->i; aj = aa->j; lu->val = aa->a;
395       } else {
396         SETERRQ(PETSC_ERR_SUP,"No mumps factorization for this matrix type");
397       }
398       if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization, get irn and jcn */
399         ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr);
400         ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr);
401         nz = 0;
402         for (i=0; i<M; i++){
403           rnz = ai[i+1] - ai[i];
404           while (rnz--) {  /* Fortran row/col index! */
405             lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++;
406           }
407         }
408       }
409     }
410     break;
411   case 3:  /* distributed assembled matrix input (size>1) */
412     if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){
413       valOnly = PETSC_FALSE;
414     } else {
415       valOnly = PETSC_TRUE; /* only update mat values, not row and col index */
416     }
417     ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr);
418     break;
419   default: SETERRQ(PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS.");
420   }
421 
422   /* analysis phase */
423   /*----------------*/
424   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){
425     lu->id.job = 1;
426 
427     lu->id.n = M;
428     switch (lu->id.ICNTL(18)){
429     case 0:  /* centralized assembled matrix input */
430       if (!lu->myid) {
431         lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn;
432         if (lu->id.ICNTL(6)>1){
433 #if defined(PETSC_USE_COMPLEX)
434           lu->id.a = (mumps_double_complex*)lu->val;
435 #else
436           lu->id.a = lu->val;
437 #endif
438         }
439       }
440       break;
441     case 3:  /* distributed assembled matrix input (size>1) */
442       lu->id.nz_loc = nnz;
443       lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn;
444       if (lu->id.ICNTL(6)>1) {
445 #if defined(PETSC_USE_COMPLEX)
446         lu->id.a_loc = (mumps_double_complex*)lu->val;
447 #else
448         lu->id.a_loc = lu->val;
449 #endif
450       }
451       /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
452       if (!lu->myid){
453         ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap.N,&lu->b_seq);CHKERRQ(ierr);
454         ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap.N,0,1,&is_iden);CHKERRQ(ierr);
455       } else {
456         ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr);
457         ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr);
458       }
459       ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr);
460       ierr = VecSetSizes(b,A->rmap.n,PETSC_DECIDE);CHKERRQ(ierr);
461       ierr = VecSetFromOptions(b);CHKERRQ(ierr);
462 
463       ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr);
464       ierr = ISDestroy(is_iden);CHKERRQ(ierr);
465       ierr = VecDestroy(b);CHKERRQ(ierr);
466       break;
467     }
468 #if defined(PETSC_USE_COMPLEX)
469     zmumps_c(&lu->id);
470 #else
471     dmumps_c(&lu->id);
472 #endif
473     if (lu->id.INFOG(1) < 0) {
474       SETERRQ1(PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1));
475     }
476   }
477 
478   /* numerical factorization phase */
479   /*-------------------------------*/
480   lu->id.job = 2;
481   if(!lu->id.ICNTL(18)) {
482     if (!lu->myid) {
483 #if defined(PETSC_USE_COMPLEX)
484       lu->id.a = (mumps_double_complex*)lu->val;
485 #else
486       lu->id.a = lu->val;
487 #endif
488     }
489   } else {
490 #if defined(PETSC_USE_COMPLEX)
491     lu->id.a_loc = (mumps_double_complex*)lu->val;
492 #else
493     lu->id.a_loc = lu->val;
494 #endif
495   }
496 #if defined(PETSC_USE_COMPLEX)
497   zmumps_c(&lu->id);
498 #else
499   dmumps_c(&lu->id);
500 #endif
501   if (lu->id.INFOG(1) < 0) {
502     if (lu->id.INFO(1) == -13) {
503       SETERRQ1(PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d megabytes\n",lu->id.INFO(2));
504     } else {
505       SETERRQ2(PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: INFO(1)=%d, INFO(2)=%d\n",lu->id.INFO(1),lu->id.INFO(2));
506     }
507   }
508 
509   if (!lu->myid && lu->id.ICNTL(16) > 0){
510     SETERRQ1(PETSC_ERR_LIB,"  lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16));
511   }
512 
513   if (lu->size > 1){
514     if ((*F)->factor == MAT_FACTOR_LU){
515       F_diag = ((Mat_MPIAIJ *)(*F)->data)->A;
516     } else {
517       F_diag = ((Mat_MPISBAIJ *)(*F)->data)->A;
518     }
519     F_diag->assembled = PETSC_TRUE;
520     if (lu->nSolve){
521       ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);
522       ierr = PetscFree(lu->id.sol_loc);CHKERRQ(ierr);
523       ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);
524     }
525   }
526   (*F)->assembled   = PETSC_TRUE;
527   lu->matstruc      = SAME_NONZERO_PATTERN;
528   lu->CleanUpMUMPS  = PETSC_TRUE;
529   lu->nSolve        = 0;
530   PetscFunctionReturn(0);
531 }
532 
533 
534 /* Note the Petsc r and c permutations are ignored */
535 #undef __FUNCT__
536 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS"
537 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
538 {
539   Mat_MUMPS      *lu = (Mat_MUMPS*)(*F)->spptr;
540 
541   PetscFunctionBegin;
542   lu->sym                 = 0;
543   lu->matstruc            = DIFFERENT_NONZERO_PATTERN;
544   PetscFunctionReturn(0);
545 }
546 
547 EXTERN_C_BEGIN
548 /*
549     The seq and mpi versions of this function are the same
550 */
551 #undef __FUNCT__
552 #define __FUNCT__ "MatGetFactor_seqaij_mumps"
553 PetscErrorCode MatGetFactor_seqaij_mumps(Mat A,MatFactorType ftype,Mat *F)
554 {
555   Mat            B;
556   PetscErrorCode ierr;
557   Mat_MUMPS      *mumps;
558 
559   PetscFunctionBegin;
560   if (ftype != MAT_FACTOR_LU) {
561     SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with MUMPS Cholesky, use SBAIJ matrix");
562   }
563   /* Create the factorization matrix */
564   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
565   ierr = MatSetSizes(B,A->rmap.n,A->cmap.n,A->rmap.N,A->cmap.N);CHKERRQ(ierr);
566   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
567   ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
568 
569   B->ops->lufactornumeric  = MatFactorNumeric_MUMPS;
570   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS;
571   B->factor                = MAT_FACTOR_LU;
572 
573   ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr);
574   mumps->CleanUpMUMPS              = PETSC_FALSE;
575   mumps->isAIJ                     = PETSC_TRUE;
576   mumps->scat_rhs                  = PETSC_NULL;
577   mumps->scat_sol                  = PETSC_NULL;
578   mumps->nSolve                    = 0;
579   mumps->MatDestroy                = A->ops->destroy;
580 
581   B->spptr                         = (void*)mumps;
582 
583   *F = B;
584   PetscFunctionReturn(0);
585 }
586 EXTERN_C_END
587 
588 EXTERN_C_BEGIN
589 #undef __FUNCT__
590 #define __FUNCT__ "MatGetFactor_mpiaij_mumps"
591 PetscErrorCode MatGetFactor_mpiaij_mumps(Mat A,MatFactorType ftype,Mat *F)
592 {
593   Mat            B;
594   PetscErrorCode ierr;
595   Mat_MUMPS      *mumps;
596 
597   PetscFunctionBegin;
598   if (ftype != MAT_FACTOR_LU) {
599     SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with MUMPS Cholesky, use SBAIJ matrix");
600   }
601   /* Create the factorization matrix */
602   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
603   ierr = MatSetSizes(B,A->rmap.n,A->cmap.n,A->rmap.N,A->cmap.N);CHKERRQ(ierr);
604   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
605   ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
606   ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
607 
608   B->ops->lufactornumeric  = MatFactorNumeric_MUMPS;
609   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS;
610   B->factor                = MAT_FACTOR_LU;
611 
612   ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr);
613   mumps->CleanUpMUMPS              = PETSC_FALSE;
614   mumps->isAIJ                     = PETSC_TRUE;
615   mumps->scat_rhs                  = PETSC_NULL;
616   mumps->scat_sol                  = PETSC_NULL;
617   mumps->nSolve                    = 0;
618   mumps->MatDestroy                = A->ops->destroy;
619 
620   B->spptr                         = (void*)mumps;
621 
622   *F = B;
623   PetscFunctionReturn(0);
624 }
625 EXTERN_C_END
626 
627 /* Note the Petsc r permutation is ignored */
628 #undef __FUNCT__
629 #define __FUNCT__ "MatCholeskyFactorSymbolic_SBAIJMUMPS"
630 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJMUMPS(Mat A,IS r,MatFactorInfo *info,Mat *F)
631 {
632   Mat_MUMPS      *lu = (Mat_MUMPS*)(*F)->spptr;
633 
634   PetscFunctionBegin;
635   lu->sym                       = 2;
636   lu->matstruc                  = DIFFERENT_NONZERO_PATTERN;
637   PetscFunctionReturn(0);
638 }
639 
640 EXTERN_C_BEGIN
641 #undef __FUNCT__
642 #define __FUNCT__ "MatGetFactor_seqsbaij_mumps"
643 PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat A,MatFactorType ftype,Mat *F)
644 {
645   Mat            B;
646   PetscErrorCode ierr;
647   Mat_MUMPS      *mumps;
648 
649   PetscFunctionBegin;
650   if (ftype != MAT_FACTOR_CHOLESKY) {
651     SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix");
652   }
653   /* Create the factorization matrix */
654   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
655   ierr = MatSetSizes(B,A->rmap.n,A->cmap.n,A->rmap.N,A->cmap.N);CHKERRQ(ierr);
656   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
657   ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr);
658   ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
659 
660   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS;
661   B->ops->choleskyfactornumeric  = MatFactorNumeric_MUMPS;
662 #if !defined(PETSC_USE_COMPLEX)
663   B->ops->getinertia             = MatGetInertia_SBAIJMUMPS;
664 #endif
665   B->factor                      = MAT_FACTOR_CHOLESKY;
666 
667   ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr);
668   mumps->CleanUpMUMPS              = PETSC_FALSE;
669   mumps->isAIJ                     = PETSC_TRUE;
670   mumps->scat_rhs                  = PETSC_NULL;
671   mumps->scat_sol                  = PETSC_NULL;
672   mumps->nSolve                    = 0;
673   mumps->MatDestroy                = A->ops->destroy;
674   B->spptr                         = (void*)mumps;
675   *F = B;
676   PetscFunctionReturn(0);
677 }
678 EXTERN_C_END
679 
680 EXTERN_C_BEGIN
681 #undef __FUNCT__
682 #define __FUNCT__ "MatGetFactor_mpisbaij_mumps"
683 PetscErrorCode MatGetFactor_mpisbaij_mumps(Mat A,MatFactorType ftype,Mat *F)
684 {
685   Mat            B;
686   PetscErrorCode ierr;
687   Mat_MUMPS      *mumps;
688 
689   PetscFunctionBegin;
690   if (ftype != MAT_FACTOR_CHOLESKY) {
691     SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix");
692   }
693   /* Create the factorization matrix */
694   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
695   ierr = MatSetSizes(B,A->rmap.n,A->cmap.n,A->rmap.N,A->cmap.N);CHKERRQ(ierr);
696   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
697   ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr);
698   ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
699 
700   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS;
701   B->ops->choleskyfactornumeric  = MatFactorNumeric_MUMPS;
702 #if !defined(PETSC_USE_COMPLEX)
703   B->ops->getinertia             = MatGetInertia_SBAIJMUMPS;
704 #endif
705   B->factor                      = MAT_FACTOR_CHOLESKY;
706 
707   ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr);
708   mumps->CleanUpMUMPS              = PETSC_FALSE;
709   mumps->isAIJ                     = PETSC_TRUE;
710   mumps->scat_rhs                  = PETSC_NULL;
711   mumps->scat_sol                  = PETSC_NULL;
712   mumps->nSolve                    = 0;
713   mumps->MatDestroy                = A->ops->destroy;
714   B->spptr                         = (void*)mumps;
715   *F = B;
716   PetscFunctionReturn(0);
717 }
718 EXTERN_C_END
719 
720 #undef __FUNCT__
721 #define __FUNCT__ "MatFactorInfo_MUMPS"
722 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) {
723   Mat_MUMPS      *lu=(Mat_MUMPS*)A->spptr;
724   PetscErrorCode ierr;
725 
726   PetscFunctionBegin;
727   /* check if matrix is mumps type */
728   if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0);
729 
730   ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr);
731   ierr = PetscViewerASCIIPrintf(viewer,"  SYM (matrix type):                  %d \n",lu->id.sym);CHKERRQ(ierr);
732   ierr = PetscViewerASCIIPrintf(viewer,"  PAR (host participation):           %d \n",lu->id.par);CHKERRQ(ierr);
733   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(1) (output for error):        %d \n",lu->id.ICNTL(1));CHKERRQ(ierr);
734   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr);
735   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(3) (output for global info):  %d \n",lu->id.ICNTL(3));CHKERRQ(ierr);
736   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(4) (level of printing):       %d \n",lu->id.ICNTL(4));CHKERRQ(ierr);
737   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(5) (input mat struct):        %d \n",lu->id.ICNTL(5));CHKERRQ(ierr);
738   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(6) (matrix prescaling):       %d \n",lu->id.ICNTL(6));CHKERRQ(ierr);
739   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(7) (matrix ordering):         %d \n",lu->id.ICNTL(7));CHKERRQ(ierr);
740   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(8) (scalling strategy):       %d \n",lu->id.ICNTL(8));CHKERRQ(ierr);
741   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(9) (A/A^T x=b is solved):     %d \n",lu->id.ICNTL(9));CHKERRQ(ierr);
742   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr);
743   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(11) (error analysis):         %d \n",lu->id.ICNTL(11));CHKERRQ(ierr);
744   if (!lu->myid && lu->id.ICNTL(11)>0) {
745     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(4) (inf norm of input mat):        %g\n",lu->id.RINFOG(4));CHKERRQ(ierr);
746     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(5) (inf norm of solution):         %g\n",lu->id.RINFOG(5));CHKERRQ(ierr);
747     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(6) (inf norm of residual):         %g\n",lu->id.RINFOG(6));CHKERRQ(ierr);
748     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",lu->id.RINFOG(7),lu->id.RINFOG(8));CHKERRQ(ierr);
749     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(9) (error estimate):               %g \n",lu->id.RINFOG(9));CHKERRQ(ierr);
750     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr);
751 
752   }
753   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(12) (efficiency control):                         %d \n",lu->id.ICNTL(12));CHKERRQ(ierr);
754   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(13) (efficiency control):                         %d \n",lu->id.ICNTL(13));CHKERRQ(ierr);
755   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr);
756   /* ICNTL(15-17) not used */
757   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(18) (input mat struct):                           %d \n",lu->id.ICNTL(18));CHKERRQ(ierr);
758   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(19) (Shur complement info):                       %d \n",lu->id.ICNTL(19));CHKERRQ(ierr);
759   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(20) (rhs sparse pattern):                         %d \n",lu->id.ICNTL(20));CHKERRQ(ierr);
760   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(21) (solution struct):                            %d \n",lu->id.ICNTL(21));CHKERRQ(ierr);
761 
762   ierr = PetscViewerASCIIPrintf(viewer,"  CNTL(1) (relative pivoting threshold):      %g \n",lu->id.CNTL(1));CHKERRQ(ierr);
763   ierr = PetscViewerASCIIPrintf(viewer,"  CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr);
764   ierr = PetscViewerASCIIPrintf(viewer,"  CNTL(3) (absolute pivoting threshold):      %g \n",lu->id.CNTL(3));CHKERRQ(ierr);
765   ierr = PetscViewerASCIIPrintf(viewer,"  CNTL(4) (value of static pivoting):         %g \n",lu->id.CNTL(4));CHKERRQ(ierr);
766 
767   /* infomation local to each processor */
768   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr);}
769   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr);
770   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
771   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr);}
772   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d]  %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr);
773   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
774   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr);}
775   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d]  %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr);
776   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
777   /*
778   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      INFO(2) (info about error or warning ): \n");CHKERRQ(ierr);}
779   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %d \n",lu->myid,lu->id.INFO(2));CHKERRQ(ierr);
780   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
781   */
782 
783   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr);}
784   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr);
785   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
786 
787   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr);}
788   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr);
789   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
790 
791   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr);}
792   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr);
793   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
794 
795   if (!lu->myid){ /* information from the host */
796     ierr = PetscViewerASCIIPrintf(viewer,"  RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr);
797     ierr = PetscViewerASCIIPrintf(viewer,"  RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr);
798     ierr = PetscViewerASCIIPrintf(viewer,"  RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr);
799 
800     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr);
801     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr);
802     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr);
803     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr);
804     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr);
805     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr);
806     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr);
807     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr);
808     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr);
809     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr);
810     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr);
811     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr);
812     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr);
813     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(16) (estimated size (in MB) of all MUMPS internal data for factorization after analysis: value on the most memory consuming processor): %d \n",lu->id.INFOG(16));CHKERRQ(ierr);
814     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(17) (estimated size of all MUMPS internal data for factorization after analysis: sum over all processors): %d \n",lu->id.INFOG(17));CHKERRQ(ierr);
815     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(18) (size of all MUMPS internal data allocated during factorization: value on the most memory consuming processor): %d \n",lu->id.INFOG(18));CHKERRQ(ierr);
816     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(19) (size of all MUMPS internal data allocated during factorization: sum over all processors): %d \n",lu->id.INFOG(19));CHKERRQ(ierr);
817      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr);
818      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(21) (size in MB of memory effectively used during factorization - value on the most memory consuming processor): %d \n",lu->id.INFOG(21));CHKERRQ(ierr);
819      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): %d \n",lu->id.INFOG(22));CHKERRQ(ierr);
820      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr);
821      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr);
822      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr);
823   }
824 
825   PetscFunctionReturn(0);
826 }
827 
828 #undef __FUNCT__
829 #define __FUNCT__ "MatView_MUMPS"
830 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer)
831 {
832   PetscErrorCode    ierr;
833   PetscTruth        iascii;
834   PetscViewerFormat format;
835 
836   PetscFunctionBegin;
837     ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
838   if (iascii) {
839     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
840     if (format == PETSC_VIEWER_ASCII_INFO){
841       ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr);
842     }
843   }
844   PetscFunctionReturn(0);
845 }
846 
847 
848 /*MC
849   MATAIJMUMPS - MATAIJMUMPS = "aijmumps" - A matrix type providing direct solvers (LU) for distributed
850   and sequential matrices via the external package MUMPS.
851 
852   If MUMPS is installed (see the manual for instructions
853   on how to declare the existence of external packages),
854   a matrix type can be constructed which invokes MUMPS solvers.
855   After calling MatCreate(...,A), simply call MatSetType(A,MATAIJMUMPS), then
856   optionally call MatSeqAIJSetPreallocation() or MatMPIAIJSetPreallocation() etc DO NOT
857   call MatCreateSeqAIJ/MPIAIJ() directly or the preallocation information will be LOST!
858 
859   If created with a single process communicator, this matrix type inherits from MATSEQAIJ.
860   Otherwise, this matrix type inherits from MATMPIAIJ.  Hence for single process communicators,
861   MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
862   for communicators controlling multiple processes.  It is recommended that you call both of
863   the above preallocation routines for simplicity.  One can also call MatConvert() for an inplace
864   conversion to or from the MATSEQAIJ or MATMPIAIJ type (depending on the communicator size)
865   without data copy AFTER the matrix values are set.
866 
867   Options Database Keys:
868 + -mat_type aijmumps - sets the matrix type to "aijmumps" during a call to MatSetFromOptions()
869 . -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric
870 . -mat_mumps_icntl_4 <0,1,2,3,4> - print level
871 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide)
872 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide)
873 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T
874 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements
875 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view
876 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide)
877 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide)
878 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide)
879 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide)
880 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold
881 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement
882 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold
883 
884   Level: beginner
885 
886 .seealso: MATSBAIJMUMPS
887 M*/
888 
889 
890 /*MC
891   MATSBAIJMUMPS - MATSBAIJMUMPS = "sbaijmumps" - A symmetric matrix type providing direct solvers (Cholesky) for
892   distributed and sequential matrices via the external package MUMPS.
893 
894   If MUMPS is installed (see the manual for instructions
895   on how to declare the existence of external packages),
896   a matrix type can be constructed which invokes MUMPS solvers.
897   After calling MatCreate(...,A), simply call MatSetType(A,MATSBAIJMUMPS), then
898   optionally call MatSeqSBAIJSetPreallocation() or MatMPISBAIJSetPreallocation() DO NOT
899   call MatCreateSeqSBAIJ/MPISBAIJ() directly or the preallocation information will be LOST!
900 
901   If created with a single process communicator, this matrix type inherits from MATSEQSBAIJ.
902   Otherwise, this matrix type inherits from MATMPISBAIJ.  Hence for single process communicators,
903   MatSeqSBAIJSetPreallocation() is supported, and similarly MatMPISBAIJSetPreallocation() is supported
904   for communicators controlling multiple processes.  It is recommended that you call both of
905   the above preallocation routines for simplicity.  One can also call MatConvert() for an inplace
906   conversion to or from the MATSEQSBAIJ or MATMPISBAIJ type (depending on the communicator size)
907   without data copy AFTER the matrix values have been set.
908 
909   Options Database Keys:
910 + -mat_type sbaijmumps - sets the matrix type to "sbaijmumps" during a call to MatSetFromOptions()
911 . -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric
912 . -mat_mumps_icntl_4 <0,...,4> - print level
913 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide)
914 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide)
915 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T
916 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements
917 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view
918 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide)
919 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide)
920 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide)
921 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide)
922 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold
923 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement
924 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold
925 
926   Level: beginner
927 
928 .seealso: MATAIJMUMPS
929 M*/
930 
931