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