xref: /petsc/src/mat/impls/aij/mpi/mumps/mumps.c (revision 3823f3589d17973688d737e11728485e97bf8777) 
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  = (mumps_double_complex*)sol_loc;
228 #else
229       lu->id.sol_loc  = 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 F,Mat A,const MatFactorInfo *info)
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 
320     /* Initialize a MUMPS instance */
321     ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid);
322     ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr);
323     lu->id.job = JOB_INIT;
324     ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr);
325     lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps);
326 
327     /* Set mumps options */
328     ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr);
329     lu->id.par=1;  /* host participates factorizaton and solve */
330     lu->id.sym=lu->sym;
331     if (lu->sym == 2){
332       ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr);
333       if (flg && icntl == 1) lu->id.sym=icntl;  /* matrix is spd */
334     }
335 #if defined(PETSC_USE_COMPLEX)
336     zmumps_c(&lu->id);
337 #else
338     dmumps_c(&lu->id);
339 #endif
340 
341     if (isSeqAIJ || isSeqSBAIJ){
342       lu->id.ICNTL(18) = 0;   /* centralized assembled matrix input */
343     } else {
344       lu->id.ICNTL(18) = 3;   /* distributed assembled matrix input */
345     }
346 
347     icntl=-1;
348     lu->id.ICNTL(4) = 1;  /* level of printing; overwrite mumps default ICNTL(4)=2 */
349     ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr);
350     if ((flg && icntl > 0) || PetscLogPrintInfo) {
351       lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */
352     } else { /* no output */
353       lu->id.ICNTL(1) = 0;  /* error message, default= 6 */
354       lu->id.ICNTL(2) = 0;  /* output stream for diagnostic printing, statistics, and warning. default=0 */
355       lu->id.ICNTL(3) = 0; /* output stream for global information, default=6 */
356     }
357     ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): column permutation and/or scaling to get a zero-free diagonal (0 to 7)","None",lu->id.ICNTL(6),&lu->id.ICNTL(6),PETSC_NULL);CHKERRQ(ierr);
358     icntl=-1;
359     ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr);
360     if (flg) {
361       if (icntl== 1){
362         SETERRQ(PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n");
363       } else {
364         lu->id.ICNTL(7) = icntl;
365       }
366     }
367     ierr = PetscOptionsInt("-mat_mumps_icntl_8","ICNTL(8): scaling strategy (-2 to 7 or 77)","None",lu->id.ICNTL(8),&lu->id.ICNTL(8),PETSC_NULL);CHKERRQ(ierr);
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): statistics related to the linear system solved (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: defines the ordering strategy with scaling constraints (0 to 3","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr);
372     ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control: with or without ScaLAPACK","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_19","ICNTL(19): Schur complement","None",lu->id.ICNTL(19),&lu->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr);
375 
376     ierr = PetscOptionsInt("-mat_mumps_icntl_22","ICNTL(22): ","None",lu->id.ICNTL(22),&lu->id.ICNTL(22),PETSC_NULL);CHKERRQ(ierr);
377     ierr = PetscOptionsInt("-mat_mumps_icntl_23","ICNTL(23): ","None",lu->id.ICNTL(23),&lu->id.ICNTL(23),PETSC_NULL);CHKERRQ(ierr);
378 
379     ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr);
380     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);
381     ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr);
382     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);
383     PetscOptionsEnd();
384   }
385 
386   /* define matrix A */
387   switch (lu->id.ICNTL(18)){
388   case 0:  /* centralized assembled matrix input (size=1) */
389     if (!lu->myid) {
390       if (isSeqAIJ){
391         Mat_SeqAIJ   *aa = (Mat_SeqAIJ*)A->data;
392         nz               = aa->nz;
393         ai = aa->i; aj = aa->j; lu->val = aa->a;
394       } else if (isSeqSBAIJ) {
395         Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data;
396         nz                  =  aa->nz;
397         ai = aa->i; aj = aa->j; lu->val = aa->a;
398       } else {
399         SETERRQ(PETSC_ERR_SUP,"No mumps factorization for this matrix type");
400       }
401       if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization, get irn and jcn */
402         ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr);
403         ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr);
404         nz = 0;
405         for (i=0; i<M; i++){
406           rnz = ai[i+1] - ai[i];
407           while (rnz--) {  /* Fortran row/col index! */
408             lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++;
409           }
410         }
411       }
412     }
413     break;
414   case 3:  /* distributed assembled matrix input (size>1) */
415     if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){
416       valOnly = PETSC_FALSE;
417     } else {
418       valOnly = PETSC_TRUE; /* only update mat values, not row and col index */
419     }
420     ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr);
421     break;
422   default: SETERRQ(PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS.");
423   }
424 
425   /* analysis phase */
426   /*----------------*/
427   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){
428     lu->id.job = 1;
429 
430     lu->id.n = M;
431     switch (lu->id.ICNTL(18)){
432     case 0:  /* centralized assembled matrix input */
433       if (!lu->myid) {
434         lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn;
435         if (lu->id.ICNTL(6)>1){
436 #if defined(PETSC_USE_COMPLEX)
437           lu->id.a = (mumps_double_complex*)lu->val;
438 #else
439           lu->id.a = lu->val;
440 #endif
441         }
442       }
443       break;
444     case 3:  /* distributed assembled matrix input (size>1) */
445       lu->id.nz_loc = nnz;
446       lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn;
447       if (lu->id.ICNTL(6)>1) {
448 #if defined(PETSC_USE_COMPLEX)
449         lu->id.a_loc = (mumps_double_complex*)lu->val;
450 #else
451         lu->id.a_loc = lu->val;
452 #endif
453       }
454       /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
455       if (!lu->myid){
456         ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr);
457         ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr);
458       } else {
459         ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr);
460         ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr);
461       }
462       ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr);
463       ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr);
464       ierr = VecSetFromOptions(b);CHKERRQ(ierr);
465 
466       ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr);
467       ierr = ISDestroy(is_iden);CHKERRQ(ierr);
468       ierr = VecDestroy(b);CHKERRQ(ierr);
469       break;
470     }
471 #if defined(PETSC_USE_COMPLEX)
472     zmumps_c(&lu->id);
473 #else
474     dmumps_c(&lu->id);
475 #endif
476     if (lu->id.INFOG(1) < 0) {
477       SETERRQ1(PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1));
478     }
479   }
480 
481   /* numerical factorization phase */
482   /*-------------------------------*/
483   lu->id.job = 2;
484   if(!lu->id.ICNTL(18)) {
485     if (!lu->myid) {
486 #if defined(PETSC_USE_COMPLEX)
487       lu->id.a = (mumps_double_complex*)lu->val;
488 #else
489       lu->id.a = lu->val;
490 #endif
491     }
492   } else {
493 #if defined(PETSC_USE_COMPLEX)
494     lu->id.a_loc = (mumps_double_complex*)lu->val;
495 #else
496     lu->id.a_loc = lu->val;
497 #endif
498   }
499 #if defined(PETSC_USE_COMPLEX)
500   zmumps_c(&lu->id);
501 #else
502   dmumps_c(&lu->id);
503 #endif
504   if (lu->id.INFOG(1) < 0) {
505     if (lu->id.INFO(1) == -13) {
506       SETERRQ1(PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d megabytes\n",lu->id.INFO(2));
507     } else {
508       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));
509     }
510   }
511 
512   if (!lu->myid && lu->id.ICNTL(16) > 0){
513     SETERRQ1(PETSC_ERR_LIB,"  lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16));
514   }
515 
516   if (lu->size > 1){
517     if ((F)->factor == MAT_FACTOR_LU){
518       F_diag = ((Mat_MPIAIJ *)(F)->data)->A;
519     } else {
520       F_diag = ((Mat_MPISBAIJ *)(F)->data)->A;
521     }
522     F_diag->assembled = PETSC_TRUE;
523     if (lu->nSolve){
524       ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);
525       ierr = PetscFree(lu->id.sol_loc);CHKERRQ(ierr);
526       ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);
527     }
528   }
529   (F)->assembled   = PETSC_TRUE;
530   lu->matstruc      = SAME_NONZERO_PATTERN;
531   lu->CleanUpMUMPS  = PETSC_TRUE;
532   lu->nSolve        = 0;
533   PetscFunctionReturn(0);
534 }
535 
536 
537 /* Note the Petsc r and c permutations are ignored */
538 #undef __FUNCT__
539 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS"
540 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
541 {
542   Mat_MUMPS      *lu = (Mat_MUMPS*)F->spptr;
543 
544   PetscFunctionBegin;
545   lu->sym                  = 0;
546   lu->matstruc             = DIFFERENT_NONZERO_PATTERN;
547   F->ops->lufactornumeric  = MatFactorNumeric_MUMPS;
548   PetscFunctionReturn(0);
549 }
550 
551 
552 /* Note the Petsc r permutation is ignored */
553 #undef __FUNCT__
554 #define __FUNCT__ "MatCholeskyFactorSymbolic_SBAIJMUMPS"
555 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJMUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info)
556 {
557   Mat_MUMPS      *lu = (Mat_MUMPS*)(F)->spptr;
558 
559   PetscFunctionBegin;
560   lu->sym                          = 2;
561   lu->matstruc                     = DIFFERENT_NONZERO_PATTERN;
562   (F)->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS;
563 #if !defined(PETSC_USE_COMPLEX)
564   (F)->ops->getinertia            = MatGetInertia_SBAIJMUMPS;
565 #endif
566   PetscFunctionReturn(0);
567 }
568 
569 #undef __FUNCT__
570 #define __FUNCT__ "MatFactorInfo_MUMPS"
571 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) {
572   Mat_MUMPS      *lu=(Mat_MUMPS*)A->spptr;
573   PetscErrorCode ierr;
574 
575   PetscFunctionBegin;
576   /* check if matrix is mumps type */
577   if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0);
578 
579   ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr);
580   ierr = PetscViewerASCIIPrintf(viewer,"  SYM (matrix type):                  %d \n",lu->id.sym);CHKERRQ(ierr);
581   ierr = PetscViewerASCIIPrintf(viewer,"  PAR (host participation):           %d \n",lu->id.par);CHKERRQ(ierr);
582   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(1) (output for error):        %d \n",lu->id.ICNTL(1));CHKERRQ(ierr);
583   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr);
584   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(3) (output for global info):  %d \n",lu->id.ICNTL(3));CHKERRQ(ierr);
585   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(4) (level of printing):       %d \n",lu->id.ICNTL(4));CHKERRQ(ierr);
586   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(5) (input mat struct):        %d \n",lu->id.ICNTL(5));CHKERRQ(ierr);
587   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(6) (matrix prescaling):       %d \n",lu->id.ICNTL(6));CHKERRQ(ierr);
588   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(7) (matrix ordering):         %d \n",lu->id.ICNTL(7));CHKERRQ(ierr);
589   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(8) (scalling strategy):       %d \n",lu->id.ICNTL(8));CHKERRQ(ierr);
590   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(9) (A/A^T x=b is solved):     %d \n",lu->id.ICNTL(9));CHKERRQ(ierr);
591   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr);
592   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(11) (error analysis):         %d \n",lu->id.ICNTL(11));CHKERRQ(ierr);
593   if (!lu->myid && lu->id.ICNTL(11)>0) {
594     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(4) (inf norm of input mat):        %g\n",lu->id.RINFOG(4));CHKERRQ(ierr);
595     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(5) (inf norm of solution):         %g\n",lu->id.RINFOG(5));CHKERRQ(ierr);
596     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(6) (inf norm of residual):         %g\n",lu->id.RINFOG(6));CHKERRQ(ierr);
597     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);
598     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(9) (error estimate):               %g \n",lu->id.RINFOG(9));CHKERRQ(ierr);
599     ierr = PetscPrintf(PETSC_COMM_SELF,"        RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr);
600 
601   }
602   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(12) (efficiency control):                         %d \n",lu->id.ICNTL(12));CHKERRQ(ierr);
603   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(13) (efficiency control):                         %d \n",lu->id.ICNTL(13));CHKERRQ(ierr);
604   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr);
605   /* ICNTL(15-17) not used */
606   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(18) (input mat struct):                           %d \n",lu->id.ICNTL(18));CHKERRQ(ierr);
607   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(19) (Shur complement info):                       %d \n",lu->id.ICNTL(19));CHKERRQ(ierr);
608   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(20) (rhs sparse pattern):                         %d \n",lu->id.ICNTL(20));CHKERRQ(ierr);
609   ierr = PetscViewerASCIIPrintf(viewer,"  ICNTL(21) (solution struct):                            %d \n",lu->id.ICNTL(21));CHKERRQ(ierr);
610 
611   ierr = PetscViewerASCIIPrintf(viewer,"  CNTL(1) (relative pivoting threshold):      %g \n",lu->id.CNTL(1));CHKERRQ(ierr);
612   ierr = PetscViewerASCIIPrintf(viewer,"  CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr);
613   ierr = PetscViewerASCIIPrintf(viewer,"  CNTL(3) (absolute pivoting threshold):      %g \n",lu->id.CNTL(3));CHKERRQ(ierr);
614   ierr = PetscViewerASCIIPrintf(viewer,"  CNTL(4) (value of static pivoting):         %g \n",lu->id.CNTL(4));CHKERRQ(ierr);
615 
616   /* infomation local to each processor */
617   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr);}
618   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr);
619   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
620   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr);}
621   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d]  %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr);
622   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
623   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr);}
624   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d]  %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr);
625   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
626   /*
627   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      INFO(2) (info about error or warning ): \n");CHKERRQ(ierr);}
628   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %d \n",lu->myid,lu->id.INFO(2));CHKERRQ(ierr);
629   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
630   */
631 
632   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);}
633   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr);
634   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
635 
636   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr);}
637   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr);
638   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
639 
640   if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, "      INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr);}
641   ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm,"             [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr);
642   ierr = PetscSynchronizedFlush(((PetscObject)A)->comm);
643 
644   if (!lu->myid){ /* information from the host */
645     ierr = PetscViewerASCIIPrintf(viewer,"  RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr);
646     ierr = PetscViewerASCIIPrintf(viewer,"  RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr);
647     ierr = PetscViewerASCIIPrintf(viewer,"  RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr);
648 
649     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr);
650     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr);
651     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr);
652     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr);
653     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr);
654     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr);
655     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr);
656     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr);
657     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr);
658     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr);
659     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr);
660     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr);
661     ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr);
662     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);
663     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);
664     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);
665     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);
666      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr);
667      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);
668      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);
669      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr);
670      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr);
671      ierr = PetscViewerASCIIPrintf(viewer,"  INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr);
672   }
673 
674   PetscFunctionReturn(0);
675 }
676 
677 #undef __FUNCT__
678 #define __FUNCT__ "MatView_MUMPS"
679 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer)
680 {
681   PetscErrorCode    ierr;
682   PetscTruth        iascii;
683   PetscViewerFormat format;
684 
685   PetscFunctionBegin;
686     ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
687   if (iascii) {
688     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
689     if (format == PETSC_VIEWER_ASCII_INFO){
690       ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr);
691     }
692   }
693   PetscFunctionReturn(0);
694 }
695 
696 
697 /*MC
698   MATAIJMUMPS - MATAIJMUMPS = "aijmumps" - A matrix type providing direct solvers (LU) for distributed
699   and sequential matrices via the external package MUMPS.
700 
701   If MUMPS is installed (see the manual for instructions
702   on how to declare the existence of external packages),
703   a matrix type can be constructed which invokes MUMPS solvers.
704   After calling MatCreate(...,A), simply call MatSetType(A,MATAIJMUMPS), then
705   optionally call MatSeqAIJSetPreallocation() or MatMPIAIJSetPreallocation() etc DO NOT
706   call MatCreateSeqAIJ/MPIAIJ() directly or the preallocation information will be LOST!
707 
708   If created with a single process communicator, this matrix type inherits from MATSEQAIJ.
709   Otherwise, this matrix type inherits from MATMPIAIJ.  Hence for single process communicators,
710   MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
711   for communicators controlling multiple processes.  It is recommended that you call both of
712   the above preallocation routines for simplicity.  One can also call MatConvert() for an inplace
713   conversion to or from the MATSEQAIJ or MATMPIAIJ type (depending on the communicator size)
714   without data copy AFTER the matrix values are set.
715 
716   Options Database Keys:
717 + -mat_type aijmumps - sets the matrix type to "aijmumps" during a call to MatSetFromOptions()
718 . -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric
719 . -mat_mumps_icntl_4 <0,1,2,3,4> - print level
720 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide)
721 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide)
722 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T
723 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements
724 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view
725 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide)
726 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide)
727 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide)
728 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide)
729 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold
730 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement
731 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold
732 
733   Level: beginner
734 
735 .seealso: MATSBAIJMUMPS
736 M*/
737 
738 
739 #undef __FUNCT__
740 #define __FUNCT__ "MatGetInfo_MUMPS"
741 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info)
742 {
743     Mat_MUMPS  *lu =(Mat_MUMPS*)A->spptr;
744 
745   PetscFunctionBegin;
746   info->block_size        = 1.0;
747   info->nz_allocated      = lu->id.INFOG(20);
748   info->nz_used           = lu->id.INFOG(20);
749   info->nz_unneeded       = 0.0;
750   info->assemblies        = 0.0;
751   info->mallocs           = 0.0;
752   info->memory            = 0.0;
753   info->fill_ratio_given  = 0;
754   info->fill_ratio_needed = 0;
755   info->factor_mallocs    = 0;
756   PetscFunctionReturn(0);
757 }
758 
759 /*MC
760   MATSBAIJMUMPS - MATSBAIJMUMPS = "sbaijmumps" - A symmetric matrix type providing direct solvers (Cholesky) for
761   distributed and sequential matrices via the external package MUMPS.
762 
763   If MUMPS is installed (see the manual for instructions
764   on how to declare the existence of external packages),
765   a matrix type can be constructed which invokes MUMPS solvers.
766   After calling MatCreate(...,A), simply call MatSetType(A,MATSBAIJMUMPS), then
767   optionally call MatSeqSBAIJSetPreallocation() or MatMPISBAIJSetPreallocation() DO NOT
768   call MatCreateSeqSBAIJ/MPISBAIJ() directly or the preallocation information will be LOST!
769 
770   If created with a single process communicator, this matrix type inherits from MATSEQSBAIJ.
771   Otherwise, this matrix type inherits from MATMPISBAIJ.  Hence for single process communicators,
772   MatSeqSBAIJSetPreallocation() is supported, and similarly MatMPISBAIJSetPreallocation() is supported
773   for communicators controlling multiple processes.  It is recommended that you call both of
774   the above preallocation routines for simplicity.  One can also call MatConvert() for an inplace
775   conversion to or from the MATSEQSBAIJ or MATMPISBAIJ type (depending on the communicator size)
776   without data copy AFTER the matrix values have been set.
777 
778   Options Database Keys:
779 + -mat_type sbaijmumps - sets the matrix type to "sbaijmumps" during a call to MatSetFromOptions()
780 . -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric
781 . -mat_mumps_icntl_4 <0,...,4> - print level
782 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide)
783 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide)
784 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T
785 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements
786 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view
787 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide)
788 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide)
789 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide)
790 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide)
791 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold
792 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement
793 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold
794 
795   Level: beginner
796 
797 .seealso: MATAIJMUMPS
798 M*/
799 
800 EXTERN_C_BEGIN
801 #undef __FUNCT__
802 #define __FUNCT__ "MatFactorGetSolverPackage_mumps"
803 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type)
804 {
805   PetscFunctionBegin;
806   *type = MAT_SOLVER_MUMPS;
807   PetscFunctionReturn(0);
808 }
809 EXTERN_C_END
810 
811 EXTERN_C_BEGIN
812 /*
813     The seq and mpi versions of this function are the same
814 */
815 #undef __FUNCT__
816 #define __FUNCT__ "MatGetFactor_seqaij_mumps"
817 PetscErrorCode MatGetFactor_seqaij_mumps(Mat A,MatFactorType ftype,Mat *F)
818 {
819   Mat            B;
820   PetscErrorCode ierr;
821   Mat_MUMPS      *mumps;
822 
823   PetscFunctionBegin;
824   if (ftype != MAT_FACTOR_LU) {
825     SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with MUMPS Cholesky, use SBAIJ matrix");
826   }
827   /* Create the factorization matrix */
828   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
829   ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
830   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
831   ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
832 
833   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS;
834   B->ops->view             = MatView_MUMPS;
835   B->ops->getinfo          = MatGetInfo_MUMPS;
836   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr);
837   B->factor                = MAT_FACTOR_LU;
838 
839   ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr);
840   mumps->CleanUpMUMPS              = PETSC_FALSE;
841   mumps->isAIJ                     = PETSC_TRUE;
842   mumps->scat_rhs                  = PETSC_NULL;
843   mumps->scat_sol                  = PETSC_NULL;
844   mumps->nSolve                    = 0;
845   mumps->MatDestroy                = B->ops->destroy;
846   B->ops->destroy                  = MatDestroy_MUMPS;
847   B->spptr                         = (void*)mumps;
848 
849   *F = B;
850   PetscFunctionReturn(0);
851 }
852 EXTERN_C_END
853 
854 EXTERN_C_BEGIN
855 #undef __FUNCT__
856 #define __FUNCT__ "MatGetFactor_mpiaij_mumps"
857 PetscErrorCode MatGetFactor_mpiaij_mumps(Mat A,MatFactorType ftype,Mat *F)
858 {
859   Mat            B;
860   PetscErrorCode ierr;
861   Mat_MUMPS      *mumps;
862 
863   PetscFunctionBegin;
864   if (ftype != MAT_FACTOR_LU) {
865     SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with MUMPS Cholesky, use SBAIJ matrix");
866   }
867   /* Create the factorization matrix */
868   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
869   ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
870   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
871   ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
872   ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
873 
874   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS;
875   B->ops->view             = MatView_MUMPS;
876   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr);
877   B->factor                = MAT_FACTOR_LU;
878 
879   ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr);
880   mumps->CleanUpMUMPS              = PETSC_FALSE;
881   mumps->isAIJ                     = PETSC_TRUE;
882   mumps->scat_rhs                  = PETSC_NULL;
883   mumps->scat_sol                  = PETSC_NULL;
884   mumps->nSolve                    = 0;
885   mumps->MatDestroy                = B->ops->destroy;
886   B->ops->destroy                  = MatDestroy_MUMPS;
887   B->spptr                         = (void*)mumps;
888 
889   *F = B;
890   PetscFunctionReturn(0);
891 }
892 EXTERN_C_END
893 
894 EXTERN_C_BEGIN
895 #undef __FUNCT__
896 #define __FUNCT__ "MatGetFactor_seqsbaij_mumps"
897 PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat A,MatFactorType ftype,Mat *F)
898 {
899   Mat            B;
900   PetscErrorCode ierr;
901   Mat_MUMPS      *mumps;
902 
903   PetscFunctionBegin;
904   if (ftype != MAT_FACTOR_CHOLESKY) {
905     SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix");
906   }
907   /* Create the factorization matrix */
908   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
909   ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
910   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
911   ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr);
912   ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
913 
914   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS;
915   B->ops->view                   = MatView_MUMPS;
916   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr);
917 
918   B->factor                      = MAT_FACTOR_CHOLESKY;
919 
920   ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr);
921   mumps->CleanUpMUMPS              = PETSC_FALSE;
922   mumps->isAIJ                     = PETSC_TRUE;
923   mumps->scat_rhs                  = PETSC_NULL;
924   mumps->scat_sol                  = PETSC_NULL;
925   mumps->nSolve                    = 0;
926   mumps->MatDestroy                = B->ops->destroy;
927   B->ops->destroy                  = MatDestroy_MUMPS;
928   B->spptr                         = (void*)mumps;
929 
930   *F = B;
931   PetscFunctionReturn(0);
932 }
933 EXTERN_C_END
934 
935 EXTERN_C_BEGIN
936 #undef __FUNCT__
937 #define __FUNCT__ "MatGetFactor_mpisbaij_mumps"
938 PetscErrorCode MatGetFactor_mpisbaij_mumps(Mat A,MatFactorType ftype,Mat *F)
939 {
940   Mat            B;
941   PetscErrorCode ierr;
942   Mat_MUMPS      *mumps;
943 
944   PetscFunctionBegin;
945   if (ftype != MAT_FACTOR_CHOLESKY) {
946     SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix");
947   }
948   /* Create the factorization matrix */
949   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
950   ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
951   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
952   ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr);
953   ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
954 
955   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS;
956   B->ops->view                   = MatView_MUMPS;
957   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr);
958   B->factor                      = MAT_FACTOR_CHOLESKY;
959 
960   ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr);
961   mumps->CleanUpMUMPS              = PETSC_FALSE;
962   mumps->isAIJ                     = PETSC_TRUE;
963   mumps->scat_rhs                  = PETSC_NULL;
964   mumps->scat_sol                  = PETSC_NULL;
965   mumps->nSolve                    = 0;
966   mumps->MatDestroy                = B->ops->destroy;
967   B->ops->destroy                  = MatDestroy_MUMPS;
968   B->spptr                         = (void*)mumps;
969 
970   *F = B;
971   PetscFunctionReturn(0);
972 }
973 EXTERN_C_END
974