Actual source code: mumps.c

petsc-master 2015-03-04
Report Typos and Errors
  2: /*
  3:     Provides an interface to the MUMPS sparse solver
  4: */

  6: #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I  "petscmat.h"  I*/
  7: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>

  9: EXTERN_C_BEGIN
 10: #if defined(PETSC_USE_COMPLEX)
 11: #if defined(PETSC_USE_REAL_SINGLE)
 12: #include <cmumps_c.h>
 13: #else
 14: #include <zmumps_c.h>
 15: #endif
 16: #else
 17: #if defined(PETSC_USE_REAL_SINGLE)
 18: #include <smumps_c.h>
 19: #else
 20: #include <dmumps_c.h>
 21: #endif
 22: #endif
 23: EXTERN_C_END
 24: #define JOB_INIT -1
 25: #define JOB_FACTSYMBOLIC 1
 26: #define JOB_FACTNUMERIC 2
 27: #define JOB_SOLVE 3
 28: #define JOB_END -2

 30: /* calls to MUMPS */
 31: #if defined(PETSC_USE_COMPLEX)
 32: #if defined(PETSC_USE_REAL_SINGLE)
 33: #define PetscMUMPS_c cmumps_c
 34: #else
 35: #define PetscMUMPS_c zmumps_c
 36: #endif
 37: #else
 38: #if defined(PETSC_USE_REAL_SINGLE)
 39: #define PetscMUMPS_c smumps_c
 40: #else
 41: #define PetscMUMPS_c dmumps_c
 42: #endif
 43: #endif

 45: /* declare MumpsScalar */
 46: #if defined(PETSC_USE_COMPLEX)
 47: #if defined(PETSC_USE_REAL_SINGLE)
 48: #define MumpsScalar mumps_complex
 49: #else
 50: #define MumpsScalar mumps_double_complex
 51: #endif
 52: #else
 53: #define MumpsScalar PetscScalar
 54: #endif

 56: /* macros s.t. indices match MUMPS documentation */
 57: #define ICNTL(I) icntl[(I)-1]
 58: #define CNTL(I) cntl[(I)-1]
 59: #define INFOG(I) infog[(I)-1]
 60: #define INFO(I) info[(I)-1]
 61: #define RINFOG(I) rinfog[(I)-1]
 62: #define RINFO(I) rinfo[(I)-1]

 64: typedef struct {
 65: #if defined(PETSC_USE_COMPLEX)
 66: #if defined(PETSC_USE_REAL_SINGLE)
 67:   CMUMPS_STRUC_C id;
 68: #else
 69:   ZMUMPS_STRUC_C id;
 70: #endif
 71: #else
 72: #if defined(PETSC_USE_REAL_SINGLE)
 73:   SMUMPS_STRUC_C id;
 74: #else
 75:   DMUMPS_STRUC_C id;
 76: #endif
 77: #endif

 79:   MatStructure matstruc;
 80:   PetscMPIInt  myid,size;
 81:   PetscInt     *irn,*jcn,nz,sym;
 82:   PetscScalar  *val;
 83:   MPI_Comm     comm_mumps;
 84:   PetscBool    isAIJ,CleanUpMUMPS;
 85:   PetscInt     ICNTL9_pre;           /* check if ICNTL(9) is changed from previous MatSolve */
 86:   VecScatter   scat_rhs, scat_sol;   /* used by MatSolve() */
 87:   Vec          b_seq,x_seq;
 88:   PetscInt     ninfo,*info;          /* display INFO */

 90:   PetscErrorCode (*Destroy)(Mat);
 91:   PetscErrorCode (*ConvertToTriples)(Mat, int, MatReuse, int*, int**, int**, PetscScalar**);
 92: } Mat_MUMPS;

 94: extern PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*);

 96: /*
 97:   MatConvertToTriples_A_B - convert Petsc matrix to triples: row[nz], col[nz], val[nz] 

 99:   input:
100:     A       - matrix in aij,baij or sbaij (bs=1) format
101:     shift   - 0: C style output triple; 1: Fortran style output triple.
102:     reuse   - MAT_INITIAL_MATRIX: spaces are allocated and values are set for the triple
103:               MAT_REUSE_MATRIX:   only the values in v array are updated
104:   output:
105:     nnz     - dim of r, c, and v (number of local nonzero entries of A)
106:     r, c, v - row and col index, matrix values (matrix triples)

108:   The returned values r, c, and sometimes v are obtained in a single PetscMalloc(). Then in MatDestroy_MUMPS() it is
109:   freed with PetscFree((mumps->irn);  This is not ideal code, the fact that v is ONLY sometimes part of mumps->irn means
110:   that the PetscMalloc() cannot easily be replaced with a PetscMalloc3(). 

112:  */

116: PetscErrorCode MatConvertToTriples_seqaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
117: {
118:   const PetscInt *ai,*aj,*ajj,M=A->rmap->n;
119:   PetscInt       nz,rnz,i,j;
121:   PetscInt       *row,*col;
122:   Mat_SeqAIJ     *aa=(Mat_SeqAIJ*)A->data;

125:   *v=aa->a;
126:   if (reuse == MAT_INITIAL_MATRIX) {
127:     nz   = aa->nz;
128:     ai   = aa->i;
129:     aj   = aa->j;
130:     *nnz = nz;
131:     PetscMalloc1(2*nz, &row);
132:     col  = row + nz;

134:     nz = 0;
135:     for (i=0; i<M; i++) {
136:       rnz = ai[i+1] - ai[i];
137:       ajj = aj + ai[i];
138:       for (j=0; j<rnz; j++) {
139:         row[nz] = i+shift; col[nz++] = ajj[j] + shift;
140:       }
141:     }
142:     *r = row; *c = col;
143:   }
144:   return(0);
145: }

149: PetscErrorCode MatConvertToTriples_seqbaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
150: {
151:   Mat_SeqBAIJ    *aa=(Mat_SeqBAIJ*)A->data;
152:   const PetscInt *ai,*aj,*ajj,bs2 = aa->bs2;
153:   PetscInt       bs,M,nz,idx=0,rnz,i,j,k,m;
155:   PetscInt       *row,*col;

158:   MatGetBlockSize(A,&bs);
159:   M = A->rmap->N/bs;
160:   *v = aa->a;
161:   if (reuse == MAT_INITIAL_MATRIX) {
162:     ai   = aa->i; aj = aa->j;
163:     nz   = bs2*aa->nz;
164:     *nnz = nz;
165:     PetscMalloc1(2*nz, &row);
166:     col  = row + nz;

168:     for (i=0; i<M; i++) {
169:       ajj = aj + ai[i];
170:       rnz = ai[i+1] - ai[i];
171:       for (k=0; k<rnz; k++) {
172:         for (j=0; j<bs; j++) {
173:           for (m=0; m<bs; m++) {
174:             row[idx]   = i*bs + m + shift;
175:             col[idx++] = bs*(ajj[k]) + j + shift;
176:           }
177:         }
178:       }
179:     }
180:     *r = row; *c = col;
181:   }
182:   return(0);
183: }

187: PetscErrorCode MatConvertToTriples_seqsbaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
188: {
189:   const PetscInt *ai, *aj,*ajj,M=A->rmap->n;
190:   PetscInt       nz,rnz,i,j;
192:   PetscInt       *row,*col;
193:   Mat_SeqSBAIJ   *aa=(Mat_SeqSBAIJ*)A->data;

196:   *v = aa->a;
197:   if (reuse == MAT_INITIAL_MATRIX) {
198:     nz   = aa->nz;
199:     ai   = aa->i;
200:     aj   = aa->j;
201:     *v   = aa->a;
202:     *nnz = nz;
203:     PetscMalloc1(2*nz, &row);
204:     col  = row + nz;

206:     nz = 0;
207:     for (i=0; i<M; i++) {
208:       rnz = ai[i+1] - ai[i];
209:       ajj = aj + ai[i];
210:       for (j=0; j<rnz; j++) {
211:         row[nz] = i+shift; col[nz++] = ajj[j] + shift;
212:       }
213:     }
214:     *r = row; *c = col;
215:   }
216:   return(0);
217: }

221: PetscErrorCode MatConvertToTriples_seqaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
222: {
223:   const PetscInt    *ai,*aj,*ajj,*adiag,M=A->rmap->n;
224:   PetscInt          nz,rnz,i,j;
225:   const PetscScalar *av,*v1;
226:   PetscScalar       *val;
227:   PetscErrorCode    ierr;
228:   PetscInt          *row,*col;
229:   Mat_SeqAIJ        *aa=(Mat_SeqAIJ*)A->data;

232:   ai   =aa->i; aj=aa->j;av=aa->a;
233:   adiag=aa->diag;
234:   if (reuse == MAT_INITIAL_MATRIX) {
235:     /* count nz in the uppper triangular part of A */
236:     nz = 0;
237:     for (i=0; i<M; i++) nz += ai[i+1] - adiag[i];
238:     *nnz = nz;

240:     PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
241:     col  = row + nz;
242:     val  = (PetscScalar*)(col + nz);

244:     nz = 0;
245:     for (i=0; i<M; i++) {
246:       rnz = ai[i+1] - adiag[i];
247:       ajj = aj + adiag[i];
248:       v1  = av + adiag[i];
249:       for (j=0; j<rnz; j++) {
250:         row[nz] = i+shift; col[nz] = ajj[j] + shift; val[nz++] = v1[j];
251:       }
252:     }
253:     *r = row; *c = col; *v = val;
254:   } else {
255:     nz = 0; val = *v;
256:     for (i=0; i <M; i++) {
257:       rnz = ai[i+1] - adiag[i];
258:       ajj = aj + adiag[i];
259:       v1  = av + adiag[i];
260:       for (j=0; j<rnz; j++) {
261:         val[nz++] = v1[j];
262:       }
263:     }
264:   }
265:   return(0);
266: }

270: PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
271: {
272:   const PetscInt    *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
273:   PetscErrorCode    ierr;
274:   PetscInt          rstart,nz,i,j,jj,irow,countA,countB;
275:   PetscInt          *row,*col;
276:   const PetscScalar *av, *bv,*v1,*v2;
277:   PetscScalar       *val;
278:   Mat_MPISBAIJ      *mat = (Mat_MPISBAIJ*)A->data;
279:   Mat_SeqSBAIJ      *aa  = (Mat_SeqSBAIJ*)(mat->A)->data;
280:   Mat_SeqBAIJ       *bb  = (Mat_SeqBAIJ*)(mat->B)->data;

283:   ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart;
284:   av=aa->a; bv=bb->a;

286:   garray = mat->garray;

288:   if (reuse == MAT_INITIAL_MATRIX) {
289:     nz   = aa->nz + bb->nz;
290:     *nnz = nz;
291:     PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
292:     col  = row + nz;
293:     val  = (PetscScalar*)(col + nz);

295:     *r = row; *c = col; *v = val;
296:   } else {
297:     row = *r; col = *c; val = *v;
298:   }

300:   jj = 0; irow = rstart;
301:   for (i=0; i<m; i++) {
302:     ajj    = aj + ai[i];                 /* ptr to the beginning of this row */
303:     countA = ai[i+1] - ai[i];
304:     countB = bi[i+1] - bi[i];
305:     bjj    = bj + bi[i];
306:     v1     = av + ai[i];
307:     v2     = bv + bi[i];

309:     /* A-part */
310:     for (j=0; j<countA; j++) {
311:       if (reuse == MAT_INITIAL_MATRIX) {
312:         row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift;
313:       }
314:       val[jj++] = v1[j];
315:     }

317:     /* B-part */
318:     for (j=0; j < countB; j++) {
319:       if (reuse == MAT_INITIAL_MATRIX) {
320:         row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift;
321:       }
322:       val[jj++] = v2[j];
323:     }
324:     irow++;
325:   }
326:   return(0);
327: }

331: PetscErrorCode MatConvertToTriples_mpiaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
332: {
333:   const PetscInt    *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
334:   PetscErrorCode    ierr;
335:   PetscInt          rstart,nz,i,j,jj,irow,countA,countB;
336:   PetscInt          *row,*col;
337:   const PetscScalar *av, *bv,*v1,*v2;
338:   PetscScalar       *val;
339:   Mat_MPIAIJ        *mat = (Mat_MPIAIJ*)A->data;
340:   Mat_SeqAIJ        *aa  = (Mat_SeqAIJ*)(mat->A)->data;
341:   Mat_SeqAIJ        *bb  = (Mat_SeqAIJ*)(mat->B)->data;

344:   ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart;
345:   av=aa->a; bv=bb->a;

347:   garray = mat->garray;

349:   if (reuse == MAT_INITIAL_MATRIX) {
350:     nz   = aa->nz + bb->nz;
351:     *nnz = nz;
352:     PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
353:     col  = row + nz;
354:     val  = (PetscScalar*)(col + nz);

356:     *r = row; *c = col; *v = val;
357:   } else {
358:     row = *r; col = *c; val = *v;
359:   }

361:   jj = 0; irow = rstart;
362:   for (i=0; i<m; i++) {
363:     ajj    = aj + ai[i];                 /* ptr to the beginning of this row */
364:     countA = ai[i+1] - ai[i];
365:     countB = bi[i+1] - bi[i];
366:     bjj    = bj + bi[i];
367:     v1     = av + ai[i];
368:     v2     = bv + bi[i];

370:     /* A-part */
371:     for (j=0; j<countA; j++) {
372:       if (reuse == MAT_INITIAL_MATRIX) {
373:         row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift;
374:       }
375:       val[jj++] = v1[j];
376:     }

378:     /* B-part */
379:     for (j=0; j < countB; j++) {
380:       if (reuse == MAT_INITIAL_MATRIX) {
381:         row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift;
382:       }
383:       val[jj++] = v2[j];
384:     }
385:     irow++;
386:   }
387:   return(0);
388: }

392: PetscErrorCode MatConvertToTriples_mpibaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
393: {
394:   Mat_MPIBAIJ       *mat    = (Mat_MPIBAIJ*)A->data;
395:   Mat_SeqBAIJ       *aa     = (Mat_SeqBAIJ*)(mat->A)->data;
396:   Mat_SeqBAIJ       *bb     = (Mat_SeqBAIJ*)(mat->B)->data;
397:   const PetscInt    *ai     = aa->i, *bi = bb->i, *aj = aa->j, *bj = bb->j,*ajj, *bjj;
398:   const PetscInt    *garray = mat->garray,mbs=mat->mbs,rstart=A->rmap->rstart;
399:   const PetscInt    bs2=mat->bs2;
400:   PetscErrorCode    ierr;
401:   PetscInt          bs,nz,i,j,k,n,jj,irow,countA,countB,idx;
402:   PetscInt          *row,*col;
403:   const PetscScalar *av=aa->a, *bv=bb->a,*v1,*v2;
404:   PetscScalar       *val;

407:   MatGetBlockSize(A,&bs);
408:   if (reuse == MAT_INITIAL_MATRIX) {
409:     nz   = bs2*(aa->nz + bb->nz);
410:     *nnz = nz;
411:     PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
412:     col  = row + nz;
413:     val  = (PetscScalar*)(col + nz);

415:     *r = row; *c = col; *v = val;
416:   } else {
417:     row = *r; col = *c; val = *v;
418:   }

420:   jj = 0; irow = rstart;
421:   for (i=0; i<mbs; i++) {
422:     countA = ai[i+1] - ai[i];
423:     countB = bi[i+1] - bi[i];
424:     ajj    = aj + ai[i];
425:     bjj    = bj + bi[i];
426:     v1     = av + bs2*ai[i];
427:     v2     = bv + bs2*bi[i];

429:     idx = 0;
430:     /* A-part */
431:     for (k=0; k<countA; k++) {
432:       for (j=0; j<bs; j++) {
433:         for (n=0; n<bs; n++) {
434:           if (reuse == MAT_INITIAL_MATRIX) {
435:             row[jj] = irow + n + shift;
436:             col[jj] = rstart + bs*ajj[k] + j + shift;
437:           }
438:           val[jj++] = v1[idx++];
439:         }
440:       }
441:     }

443:     idx = 0;
444:     /* B-part */
445:     for (k=0; k<countB; k++) {
446:       for (j=0; j<bs; j++) {
447:         for (n=0; n<bs; n++) {
448:           if (reuse == MAT_INITIAL_MATRIX) {
449:             row[jj] = irow + n + shift;
450:             col[jj] = bs*garray[bjj[k]] + j + shift;
451:           }
452:           val[jj++] = v2[idx++];
453:         }
454:       }
455:     }
456:     irow += bs;
457:   }
458:   return(0);
459: }

463: PetscErrorCode MatConvertToTriples_mpiaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
464: {
465:   const PetscInt    *ai, *aj,*adiag, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
466:   PetscErrorCode    ierr;
467:   PetscInt          rstart,nz,nza,nzb,i,j,jj,irow,countA,countB;
468:   PetscInt          *row,*col;
469:   const PetscScalar *av, *bv,*v1,*v2;
470:   PetscScalar       *val;
471:   Mat_MPIAIJ        *mat =  (Mat_MPIAIJ*)A->data;
472:   Mat_SeqAIJ        *aa  =(Mat_SeqAIJ*)(mat->A)->data;
473:   Mat_SeqAIJ        *bb  =(Mat_SeqAIJ*)(mat->B)->data;

476:   ai=aa->i; aj=aa->j; adiag=aa->diag;
477:   bi=bb->i; bj=bb->j; garray = mat->garray;
478:   av=aa->a; bv=bb->a;

480:   rstart = A->rmap->rstart;

482:   if (reuse == MAT_INITIAL_MATRIX) {
483:     nza = 0;    /* num of upper triangular entries in mat->A, including diagonals */
484:     nzb = 0;    /* num of upper triangular entries in mat->B */
485:     for (i=0; i<m; i++) {
486:       nza   += (ai[i+1] - adiag[i]);
487:       countB = bi[i+1] - bi[i];
488:       bjj    = bj + bi[i];
489:       for (j=0; j<countB; j++) {
490:         if (garray[bjj[j]] > rstart) nzb++;
491:       }
492:     }

494:     nz   = nza + nzb; /* total nz of upper triangular part of mat */
495:     *nnz = nz;
496:     PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
497:     col  = row + nz;
498:     val  = (PetscScalar*)(col + nz);

500:     *r = row; *c = col; *v = val;
501:   } else {
502:     row = *r; col = *c; val = *v;
503:   }

505:   jj = 0; irow = rstart;
506:   for (i=0; i<m; i++) {
507:     ajj    = aj + adiag[i];                 /* ptr to the beginning of the diagonal of this row */
508:     v1     = av + adiag[i];
509:     countA = ai[i+1] - adiag[i];
510:     countB = bi[i+1] - bi[i];
511:     bjj    = bj + bi[i];
512:     v2     = bv + bi[i];

514:     /* A-part */
515:     for (j=0; j<countA; j++) {
516:       if (reuse == MAT_INITIAL_MATRIX) {
517:         row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift;
518:       }
519:       val[jj++] = v1[j];
520:     }

522:     /* B-part */
523:     for (j=0; j < countB; j++) {
524:       if (garray[bjj[j]] > rstart) {
525:         if (reuse == MAT_INITIAL_MATRIX) {
526:           row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift;
527:         }
528:         val[jj++] = v2[j];
529:       }
530:     }
531:     irow++;
532:   }
533:   return(0);
534: }

538: PetscErrorCode MatGetDiagonal_MUMPS(Mat A,Vec v)
539: {
541:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: MUMPS factor");
542:   return(0);
543: }

547: PetscErrorCode MatDestroy_MUMPS(Mat A)
548: {
549:   Mat_MUMPS      *mumps=(Mat_MUMPS*)A->spptr;

553:   if (mumps->CleanUpMUMPS) {
554:     /* Terminate instance, deallocate memories */
555:     PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);
556:     VecScatterDestroy(&mumps->scat_rhs);
557:     VecScatterDestroy(&mumps->scat_sol);
558:     VecDestroy(&mumps->b_seq);
559:     VecDestroy(&mumps->x_seq);
560:     PetscFree(mumps->id.perm_in);
561:     PetscFree(mumps->irn);
562:     PetscFree(mumps->info);

564:     mumps->id.job = JOB_END;
565:     PetscMUMPS_c(&mumps->id);
566:     MPI_Comm_free(&(mumps->comm_mumps));
567:   }
568:   if (mumps->Destroy) {
569:     (mumps->Destroy)(A);
570:   }
571:   PetscFree(A->spptr);

573:   /* clear composed functions */
574:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
575:   PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetIcntl_C",NULL);
576:   PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetIcntl_C",NULL);
577:   PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetCntl_C",NULL);
578:   PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetCntl_C",NULL);

580:   PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfo_C",NULL);
581:   PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfog_C",NULL);
582:   PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfo_C",NULL);
583:   PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfog_C",NULL);
584:   return(0);
585: }

589: PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x)
590: {
591:   Mat_MUMPS        *mumps=(Mat_MUMPS*)A->spptr;
592:   PetscScalar      *array;
593:   Vec              b_seq;
594:   IS               is_iden,is_petsc;
595:   PetscErrorCode   ierr;
596:   PetscInt         i;
597:   static PetscBool cite1 = PETSC_FALSE,cite2 = PETSC_FALSE;

600:   PetscCitationsRegister("@article{MUMPS01,\n  author = {P.~R. Amestoy and I.~S. Duff and J.-Y. L'Excellent and J. Koster},\n  title = {A fully asynchronous multifrontal solver using distributed dynamic scheduling},\n  journal = {SIAM Journal on Matrix Analysis and Applications},\n  volume = {23},\n  number = {1},\n  pages = {15--41},\n  year = {2001}\n}\n",&cite1);
601:   PetscCitationsRegister("@article{MUMPS02,\n  author = {P.~R. Amestoy and A. Guermouche and J.-Y. L'Excellent and S. Pralet},\n  title = {Hybrid scheduling for the parallel solution of linear systems},\n  journal = {Parallel Computing},\n  volume = {32},\n  number = {2},\n  pages = {136--156},\n  year = {2006}\n}\n",&cite2);
602:   mumps->id.nrhs = 1;
603:   b_seq          = mumps->b_seq;
604:   if (mumps->size > 1) {
605:     /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */
606:     VecScatterBegin(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);
607:     VecScatterEnd(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);
608:     if (!mumps->myid) {VecGetArray(b_seq,&array);}
609:   } else {  /* size == 1 */
610:     VecCopy(b,x);
611:     VecGetArray(x,&array);
612:   }
613:   if (!mumps->myid) { /* define rhs on the host */
614:     mumps->id.nrhs = 1;
615:     mumps->id.rhs = (MumpsScalar*)array;
616:   }

618:   /* solve phase */
619:   /*-------------*/
620:   mumps->id.job = JOB_SOLVE;
621:   PetscMUMPS_c(&mumps->id);
622:   if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",mumps->id.INFOG(1));

624:   if (mumps->size > 1) { /* convert mumps distributed solution to petsc mpi x */
625:     if (mumps->scat_sol && mumps->ICNTL9_pre != mumps->id.ICNTL(9)) {
626:       /* when id.ICNTL(9) changes, the contents of lsol_loc may change (not its size, lsol_loc), recreates scat_sol */
627:       VecScatterDestroy(&mumps->scat_sol);
628:     }
629:     if (!mumps->scat_sol) { /* create scatter scat_sol */
630:       ISCreateStride(PETSC_COMM_SELF,mumps->id.lsol_loc,0,1,&is_iden); /* from */
631:       for (i=0; i<mumps->id.lsol_loc; i++) {
632:         mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */
633:       }
634:       ISCreateGeneral(PETSC_COMM_SELF,mumps->id.lsol_loc,mumps->id.isol_loc,PETSC_COPY_VALUES,&is_petsc);  /* to */
635:       VecScatterCreate(mumps->x_seq,is_iden,x,is_petsc,&mumps->scat_sol);
636:       ISDestroy(&is_iden);
637:       ISDestroy(&is_petsc);

639:       mumps->ICNTL9_pre = mumps->id.ICNTL(9); /* save current value of id.ICNTL(9) */
640:     }

642:     VecScatterBegin(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
643:     VecScatterEnd(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
644:   }
645:   return(0);
646: }

650: PetscErrorCode MatSolveTranspose_MUMPS(Mat A,Vec b,Vec x)
651: {
652:   Mat_MUMPS      *mumps=(Mat_MUMPS*)A->spptr;

656:   mumps->id.ICNTL(9) = 0;
657:   MatSolve_MUMPS(A,b,x);
658:   mumps->id.ICNTL(9) = 1;
659:   return(0);
660: }

664: PetscErrorCode MatMatSolve_MUMPS(Mat A,Mat B,Mat X)
665: {
667:   PetscBool      flg;
668:   Mat_MUMPS      *mumps=(Mat_MUMPS*)A->spptr;
669:   PetscInt       i,nrhs,M;
670:   PetscScalar    *array,*bray;

673:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
674:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
675:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
676:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
677:   if (B->rmap->n != X->rmap->n) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONG,"Matrix B and X must have same row distribution");

679:   MatGetSize(B,&M,&nrhs);
680:   mumps->id.nrhs = nrhs;
681:   mumps->id.lrhs = M;

683:   if (mumps->size == 1) {
684:     /* copy B to X */
685:     MatDenseGetArray(B,&bray);
686:     MatDenseGetArray(X,&array);
687:     for (i=0; i<M*nrhs; i++) array[i] = bray[i];
688:     MatDenseRestoreArray(B,&bray);
689:     mumps->id.rhs = (MumpsScalar*)array;

691:     /* solve phase */
692:     /*-------------*/
693:     mumps->id.job = JOB_SOLVE;
694:     PetscMUMPS_c(&mumps->id);
695:     if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",mumps->id.INFOG(1));
696:     MatDenseRestoreArray(X,&array);
697:   } else {  /*--------- parallel case --------*/
698:     PetscInt       lsol_loc,nlsol_loc,*isol_loc,*idx,*iidx,*idxx,*isol_loc_save;
699:     MumpsScalar    *sol_loc,*sol_loc_save;
700:     IS             is_to,is_from;
701:     PetscInt       k,proc,j,m;
702:     const PetscInt *rstart;
703:     Vec            v_mpi,b_seq,x_seq;
704:     VecScatter     scat_rhs,scat_sol;

706:     /* create x_seq to hold local solution */
707:     isol_loc_save = mumps->id.isol_loc; /* save it for MatSovle() */
708:     sol_loc_save  = mumps->id.sol_loc;

710:     lsol_loc  = mumps->id.INFO(23);
711:     nlsol_loc = nrhs*lsol_loc;     /* length of sol_loc */
712:     PetscMalloc2(nlsol_loc,&sol_loc,nlsol_loc,&isol_loc);
713:     mumps->id.sol_loc = (MumpsScalar*)sol_loc;
714:     mumps->id.isol_loc = isol_loc;

716:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,nlsol_loc,(PetscScalar*)sol_loc,&x_seq);

718:     /* copy rhs matrix B into vector v_mpi */
719:     MatGetLocalSize(B,&m,NULL);
720:     MatDenseGetArray(B,&bray);
721:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)B),1,nrhs*m,nrhs*M,(const PetscScalar*)bray,&v_mpi);
722:     MatDenseRestoreArray(B,&bray);

724:     /* scatter v_mpi to b_seq because MUMPS only supports centralized rhs */
725:     /* idx: maps from k-th index of v_mpi to (i,j)-th global entry of B;
726:       iidx: inverse of idx, will be used by scattering xx_seq -> X       */
727:     PetscMalloc2(nrhs*M,&idx,nrhs*M,&iidx);
728:     MatGetOwnershipRanges(B,&rstart);
729:     k = 0;
730:     for (proc=0; proc<mumps->size; proc++){
731:       for (j=0; j<nrhs; j++){
732:         for (i=rstart[proc]; i<rstart[proc+1]; i++){
733:           iidx[j*M + i] = k;
734:           idx[k++]      = j*M + i;
735:         }
736:       }
737:     }

739:     if (!mumps->myid) {
740:       VecCreateSeq(PETSC_COMM_SELF,nrhs*M,&b_seq);
741:       ISCreateGeneral(PETSC_COMM_SELF,nrhs*M,idx,PETSC_COPY_VALUES,&is_to);
742:       ISCreateStride(PETSC_COMM_SELF,nrhs*M,0,1,&is_from);
743:     } else {
744:       VecCreateSeq(PETSC_COMM_SELF,0,&b_seq);
745:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_to);
746:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_from);
747:     }
748:     VecScatterCreate(v_mpi,is_from,b_seq,is_to,&scat_rhs);
749:     VecScatterBegin(scat_rhs,v_mpi,b_seq,INSERT_VALUES,SCATTER_FORWARD);
750:     ISDestroy(&is_to);
751:     ISDestroy(&is_from);
752:     VecScatterEnd(scat_rhs,v_mpi,b_seq,INSERT_VALUES,SCATTER_FORWARD);

754:     if (!mumps->myid) { /* define rhs on the host */
755:       VecGetArray(b_seq,&bray);
756:       mumps->id.rhs = (MumpsScalar*)bray;
757:       VecRestoreArray(b_seq,&bray);
758:     }

760:     /* solve phase */
761:     /*-------------*/
762:     mumps->id.job = JOB_SOLVE;
763:     PetscMUMPS_c(&mumps->id);
764:     if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",mumps->id.INFOG(1));

766:     /* scatter mumps distributed solution to petsc vector v_mpi, which shares local arrays with solution matrix X */
767:     MatDenseGetArray(X,&array);
768:     VecPlaceArray(v_mpi,array);
769: 
770:     /* create scatter scat_sol */
771:     PetscMalloc1(nlsol_loc,&idxx);
772:     ISCreateStride(PETSC_COMM_SELF,nlsol_loc,0,1,&is_from);
773:     for (i=0; i<lsol_loc; i++) {
774:       isol_loc[i] -= 1; /* change Fortran style to C style */
775:       idxx[i] = iidx[isol_loc[i]];
776:       for (j=1; j<nrhs; j++){
777:         idxx[j*lsol_loc+i] = iidx[isol_loc[i]+j*M];
778:       }
779:     }
780:     ISCreateGeneral(PETSC_COMM_SELF,nlsol_loc,idxx,PETSC_COPY_VALUES,&is_to);
781:     VecScatterCreate(x_seq,is_from,v_mpi,is_to,&scat_sol);
782:     VecScatterBegin(scat_sol,x_seq,v_mpi,INSERT_VALUES,SCATTER_FORWARD);
783:     ISDestroy(&is_from);
784:     ISDestroy(&is_to);
785:     VecScatterEnd(scat_sol,x_seq,v_mpi,INSERT_VALUES,SCATTER_FORWARD);
786:     MatDenseRestoreArray(X,&array);

788:     /* free spaces */
789:     mumps->id.sol_loc = sol_loc_save;
790:     mumps->id.isol_loc = isol_loc_save;

792:     PetscFree2(sol_loc,isol_loc);
793:     PetscFree2(idx,iidx);
794:     PetscFree(idxx);
795:     VecDestroy(&x_seq);
796:     VecDestroy(&v_mpi);
797:     VecDestroy(&b_seq);
798:     VecScatterDestroy(&scat_rhs);
799:     VecScatterDestroy(&scat_sol);
800:   }
801:   return(0);
802: }

804: #if !defined(PETSC_USE_COMPLEX)
805: /*
806:   input:
807:    F:        numeric factor
808:   output:
809:    nneg:     total number of negative pivots
810:    nzero:    0
811:    npos:     (global dimension of F) - nneg
812: */

816: PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos)
817: {
818:   Mat_MUMPS      *mumps =(Mat_MUMPS*)F->spptr;
820:   PetscMPIInt    size;

823:   MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);
824:   /* 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 */
825:   if (size > 1 && mumps->id.ICNTL(13) != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"ICNTL(13)=%d. -mat_mumps_icntl_13 must be set as 1 for correct global matrix inertia\n",mumps->id.INFOG(13));

827:   if (nneg) *nneg = mumps->id.INFOG(12);
828:   if (nzero || npos) {
829:     if (mumps->id.ICNTL(24) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"-mat_mumps_icntl_24 must be set as 1 for null pivot row detection");
830:     if (nzero) *nzero = mumps->id.INFOG(28);
831:     if (npos) *npos   = F->rmap->N - (mumps->id.INFOG(12) + mumps->id.INFOG(28));
832:   }
833:   return(0);
834: }
835: #endif /* !defined(PETSC_USE_COMPLEX) */

839: PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info)
840: {
841:   Mat_MUMPS      *mumps =(Mat_MUMPS*)(F)->spptr;
843:   Mat            F_diag;
844:   PetscBool      isMPIAIJ;

847:   (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);

849:   /* numerical factorization phase */
850:   /*-------------------------------*/
851:   mumps->id.job = JOB_FACTNUMERIC;
852:   if (!mumps->id.ICNTL(18)) { /* A is centralized */
853:     if (!mumps->myid) {
854:       mumps->id.a = (MumpsScalar*)mumps->val;
855:     }
856:   } else {
857:     mumps->id.a_loc = (MumpsScalar*)mumps->val;
858:   }
859:   PetscMUMPS_c(&mumps->id);
860:   if (mumps->id.INFOG(1) < 0) {
861:     if (mumps->id.INFO(1) == -13) {
862:       if (mumps->id.INFO(2) < 0) {
863:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d megabytes\n",-mumps->id.INFO(2));
864:       } else {
865:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d bytes\n",mumps->id.INFO(2));
866:       }
867:     } else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: INFO(1)=%d, INFO(2)=%d\n",mumps->id.INFO(1),mumps->id.INFO(2));
868:   }
869:   if (!mumps->myid && mumps->id.ICNTL(16) > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"  mumps->id.ICNTL(16):=%d\n",mumps->id.INFOG(16));
870: 
871:   (F)->assembled      = PETSC_TRUE;
872:   mumps->matstruc     = SAME_NONZERO_PATTERN;
873:   mumps->CleanUpMUMPS = PETSC_TRUE;

875:   if (mumps->size > 1) {
876:     PetscInt    lsol_loc;
877:     PetscScalar *sol_loc;

879:     PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
880:     if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A;
881:     else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A;
882:     F_diag->assembled = PETSC_TRUE;

884:     /* distributed solution; Create x_seq=sol_loc for repeated use */
885:     if (mumps->x_seq) {
886:       VecScatterDestroy(&mumps->scat_sol);
887:       PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);
888:       VecDestroy(&mumps->x_seq);
889:     }
890:     lsol_loc = mumps->id.INFO(23); /* length of sol_loc */
891:     PetscMalloc2(lsol_loc,&sol_loc,lsol_loc,&mumps->id.isol_loc);
892:     mumps->id.lsol_loc = lsol_loc;
893:     mumps->id.sol_loc = (MumpsScalar*)sol_loc;
894:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);
895:   }
896:   return(0);
897: }

899: /* Sets MUMPS options from the options database */
902: PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A)
903: {
904:   Mat_MUMPS      *mumps = (Mat_MUMPS*)F->spptr;
906:   PetscInt       icntl,info[40],i,ninfo=40;
907:   PetscBool      flg;

910:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");
911:   PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);
912:   if (flg) mumps->id.ICNTL(1) = icntl;
913:   PetscOptionsInt("-mat_mumps_icntl_2","ICNTL(2): output stream for diagnostic printing, statistics, and warning","None",mumps->id.ICNTL(2),&icntl,&flg);
914:   if (flg) mumps->id.ICNTL(2) = icntl;
915:   PetscOptionsInt("-mat_mumps_icntl_3","ICNTL(3): output stream for global information, collected on the host","None",mumps->id.ICNTL(3),&icntl,&flg);
916:   if (flg) mumps->id.ICNTL(3) = icntl;

918:   PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);
919:   if (flg) mumps->id.ICNTL(4) = icntl;
920:   if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */

922:   PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): permutes to a zero-free diagonal and/or scale the matrix (0 to 7)","None",mumps->id.ICNTL(6),&icntl,&flg);
923:   if (flg) mumps->id.ICNTL(6) = icntl;

925:   PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): computes a symmetric permutation in sequential analysis (0 to 7). 3=Scotch, 4=PORD, 5=Metis","None",mumps->id.ICNTL(7),&icntl,&flg);
926:   if (flg) {
927:     if (icntl== 1 && mumps->size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n");
928:     else mumps->id.ICNTL(7) = icntl;
929:   }

931:   PetscOptionsInt("-mat_mumps_icntl_8","ICNTL(8): scaling strategy (-2 to 8 or 77)","None",mumps->id.ICNTL(8),&mumps->id.ICNTL(8),NULL);
932:   /* PetscOptionsInt("-mat_mumps_icntl_9","ICNTL(9): computes the solution using A or A^T","None",mumps->id.ICNTL(9),&mumps->id.ICNTL(9),NULL); handled by MatSolveTranspose_MUMPS() */
933:   PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);
934:   PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): statistics related to an error analysis (via -ksp_view)","None",mumps->id.ICNTL(11),&mumps->id.ICNTL(11),NULL);
935:   PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): an ordering strategy for symmetric matrices (0 to 3)","None",mumps->id.ICNTL(12),&mumps->id.ICNTL(12),NULL);
936:   PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): parallelism of the root node (enable ScaLAPACK) and its splitting","None",mumps->id.ICNTL(13),&mumps->id.ICNTL(13),NULL);
937:   PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): percentage increase in the estimated working space","None",mumps->id.ICNTL(14),&mumps->id.ICNTL(14),NULL);
938:   PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): computes the Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);
939:   /* PetscOptionsInt("-mat_mumps_icntl_20","ICNTL(20): the format (dense or sparse) of the right-hand sides","None",mumps->id.ICNTL(20),&mumps->id.ICNTL(20),NULL); -- sparse rhs is not supported in PETSc API */
940:   /* PetscOptionsInt("-mat_mumps_icntl_21","ICNTL(21): the distribution (centralized or distributed) of the solution vectors","None",mumps->id.ICNTL(21),&mumps->id.ICNTL(21),NULL); we only use distributed solution vector */

942:   PetscOptionsInt("-mat_mumps_icntl_22","ICNTL(22): in-core/out-of-core factorization and solve (0 or 1)","None",mumps->id.ICNTL(22),&mumps->id.ICNTL(22),NULL);
943:   PetscOptionsInt("-mat_mumps_icntl_23","ICNTL(23): max size of the working memory (MB) that can allocate per processor","None",mumps->id.ICNTL(23),&mumps->id.ICNTL(23),NULL);
944:   PetscOptionsInt("-mat_mumps_icntl_24","ICNTL(24): detection of null pivot rows (0 or 1)","None",mumps->id.ICNTL(24),&mumps->id.ICNTL(24),NULL);
945:   if (mumps->id.ICNTL(24)) {
946:     mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */
947:   }

949:   PetscOptionsInt("-mat_mumps_icntl_25","ICNTL(25): compute a solution of a deficient matrix and a null space basis","None",mumps->id.ICNTL(25),&mumps->id.ICNTL(25),NULL);
950:   PetscOptionsInt("-mat_mumps_icntl_26","ICNTL(26): drives the solution phase if a Schur complement matrix","None",mumps->id.ICNTL(26),&mumps->id.ICNTL(26),NULL);
951:   PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): the blocking size for multiple right-hand sides","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);
952:   PetscOptionsInt("-mat_mumps_icntl_28","ICNTL(28): use 1 for sequential analysis and ictnl(7) ordering, or 2 for parallel analysis and ictnl(29) ordering","None",mumps->id.ICNTL(28),&mumps->id.ICNTL(28),NULL);
953:   PetscOptionsInt("-mat_mumps_icntl_29","ICNTL(29): parallel ordering 1 = ptscotch, 2 = parmetis","None",mumps->id.ICNTL(29),&mumps->id.ICNTL(29),NULL);
954:   PetscOptionsInt("-mat_mumps_icntl_30","ICNTL(30): compute user-specified set of entries in inv(A)","None",mumps->id.ICNTL(30),&mumps->id.ICNTL(30),NULL);
955:   PetscOptionsInt("-mat_mumps_icntl_31","ICNTL(31): indicates which factors may be discarded during factorization","None",mumps->id.ICNTL(31),&mumps->id.ICNTL(31),NULL);
956:   /* PetscOptionsInt("-mat_mumps_icntl_32","ICNTL(32): performs the forward elemination of the right-hand sides during factorization","None",mumps->id.ICNTL(32),&mumps->id.ICNTL(32),NULL);  -- not supported by PETSc API */
957:   PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);

959:   PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);
960:   PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);
961:   PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);
962:   PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);
963:   PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);

965:   PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL);

967:   PetscOptionsIntArray("-mat_mumps_view_info","request INFO local to each processor","",info,&ninfo,NULL);
968:   if (ninfo) {
969:     if (ninfo > 40) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"number of INFO %d must <= 40\n",ninfo);
970:     PetscMalloc1(ninfo,&mumps->info);
971:     mumps->ninfo = ninfo;
972:     for (i=0; i<ninfo; i++) {
973:       if (info[i] < 0 || info[i]>40) {
974:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index of INFO %d must between 1 and 40\n",ninfo);
975:       } else {
976:         mumps->info[i] = info[i];
977:       }
978:     }
979:   }

981:   PetscOptionsEnd();
982:   return(0);
983: }

987: PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps)
988: {

992:   MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid);
993:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);
994:   MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));

996:   mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps);

998:   mumps->id.job = JOB_INIT;
999:   mumps->id.par = 1;  /* host participates factorizaton and solve */
1000:   mumps->id.sym = mumps->sym;
1001:   PetscMUMPS_c(&mumps->id);

1003:   mumps->CleanUpMUMPS = PETSC_FALSE;
1004:   mumps->scat_rhs     = NULL;
1005:   mumps->scat_sol     = NULL;

1007:   /* set PETSc-MUMPS default options - override MUMPS default */
1008:   mumps->id.ICNTL(3) = 0;
1009:   mumps->id.ICNTL(4) = 0;
1010:   if (mumps->size == 1) {
1011:     mumps->id.ICNTL(18) = 0;   /* centralized assembled matrix input */
1012:   } else {
1013:     mumps->id.ICNTL(18) = 3;   /* distributed assembled matrix input */
1014:     mumps->id.ICNTL(20) = 0;   /* rhs is in dense format */
1015:     mumps->id.ICNTL(21) = 1;   /* distributed solution */
1016:   }
1017:   return(0);
1018: }

1020: /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */
1023: PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
1024: {
1025:   Mat_MUMPS      *mumps = (Mat_MUMPS*)F->spptr;
1027:   Vec            b;
1028:   IS             is_iden;
1029:   const PetscInt M = A->rmap->N;

1032:   mumps->matstruc = DIFFERENT_NONZERO_PATTERN;

1034:   /* Set MUMPS options from the options database */
1035:   PetscSetMUMPSFromOptions(F,A);

1037:   (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);

1039:   /* analysis phase */
1040:   /*----------------*/
1041:   mumps->id.job = JOB_FACTSYMBOLIC;
1042:   mumps->id.n   = M;
1043:   switch (mumps->id.ICNTL(18)) {
1044:   case 0:  /* centralized assembled matrix input */
1045:     if (!mumps->myid) {
1046:       mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
1047:       if (mumps->id.ICNTL(6)>1) {
1048:         mumps->id.a = (MumpsScalar*)mumps->val;
1049:       }
1050:       if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */
1051:         /*
1052:         PetscBool      flag;
1053:         ISEqual(r,c,&flag);
1054:         if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm");
1055:         ISView(r,PETSC_VIEWER_STDOUT_SELF);
1056:          */
1057:         if (!mumps->myid) {
1058:           const PetscInt *idx;
1059:           PetscInt       i,*perm_in;

1061:           PetscMalloc1(M,&perm_in);
1062:           ISGetIndices(r,&idx);

1064:           mumps->id.perm_in = perm_in;
1065:           for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */
1066:           ISRestoreIndices(r,&idx);
1067:         }
1068:       }
1069:     }
1070:     break;
1071:   case 3:  /* distributed assembled matrix input (size>1) */
1072:     mumps->id.nz_loc = mumps->nz;
1073:     mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
1074:     if (mumps->id.ICNTL(6)>1) {
1075:       mumps->id.a_loc = (MumpsScalar*)mumps->val;
1076:     }
1077:     /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
1078:     if (!mumps->myid) {
1079:       VecCreateSeq(PETSC_COMM_SELF,A->rmap->N,&mumps->b_seq);
1080:       ISCreateStride(PETSC_COMM_SELF,A->rmap->N,0,1,&is_iden);
1081:     } else {
1082:       VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
1083:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
1084:     }
1085:     MatCreateVecs(A,NULL,&b);
1086:     VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
1087:     ISDestroy(&is_iden);
1088:     VecDestroy(&b);
1089:     break;
1090:   }
1091:   PetscMUMPS_c(&mumps->id);
1092:   if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1));

1094:   F->ops->lufactornumeric = MatFactorNumeric_MUMPS;
1095:   F->ops->solve           = MatSolve_MUMPS;
1096:   F->ops->solvetranspose  = MatSolveTranspose_MUMPS;
1097:   F->ops->matsolve        = MatMatSolve_MUMPS;
1098:   return(0);
1099: }

1101: /* Note the Petsc r and c permutations are ignored */
1104: PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
1105: {
1106:   Mat_MUMPS      *mumps = (Mat_MUMPS*)F->spptr;
1108:   Vec            b;
1109:   IS             is_iden;
1110:   const PetscInt M = A->rmap->N;

1113:   mumps->matstruc = DIFFERENT_NONZERO_PATTERN;

1115:   /* Set MUMPS options from the options database */
1116:   PetscSetMUMPSFromOptions(F,A);

1118:   (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);

1120:   /* analysis phase */
1121:   /*----------------*/
1122:   mumps->id.job = JOB_FACTSYMBOLIC;
1123:   mumps->id.n   = M;
1124:   switch (mumps->id.ICNTL(18)) {
1125:   case 0:  /* centralized assembled matrix input */
1126:     if (!mumps->myid) {
1127:       mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
1128:       if (mumps->id.ICNTL(6)>1) {
1129:         mumps->id.a = (MumpsScalar*)mumps->val;
1130:       }
1131:     }
1132:     break;
1133:   case 3:  /* distributed assembled matrix input (size>1) */
1134:     mumps->id.nz_loc = mumps->nz;
1135:     mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
1136:     if (mumps->id.ICNTL(6)>1) {
1137:       mumps->id.a_loc = (MumpsScalar*)mumps->val;
1138:     }
1139:     /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
1140:     if (!mumps->myid) {
1141:       VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);
1142:       ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
1143:     } else {
1144:       VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
1145:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
1146:     }
1147:     MatCreateVecs(A,NULL,&b);
1148:     VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
1149:     ISDestroy(&is_iden);
1150:     VecDestroy(&b);
1151:     break;
1152:   }
1153:   PetscMUMPS_c(&mumps->id);
1154:   if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1));

1156:   F->ops->lufactornumeric = MatFactorNumeric_MUMPS;
1157:   F->ops->solve           = MatSolve_MUMPS;
1158:   F->ops->solvetranspose  = MatSolveTranspose_MUMPS;
1159:   return(0);
1160: }

1162: /* Note the Petsc r permutation and factor info are ignored */
1165: PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info)
1166: {
1167:   Mat_MUMPS      *mumps = (Mat_MUMPS*)F->spptr;
1169:   Vec            b;
1170:   IS             is_iden;
1171:   const PetscInt M = A->rmap->N;

1174:   mumps->matstruc = DIFFERENT_NONZERO_PATTERN;

1176:   /* Set MUMPS options from the options database */
1177:   PetscSetMUMPSFromOptions(F,A);

1179:   (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);

1181:   /* analysis phase */
1182:   /*----------------*/
1183:   mumps->id.job = JOB_FACTSYMBOLIC;
1184:   mumps->id.n   = M;
1185:   switch (mumps->id.ICNTL(18)) {
1186:   case 0:  /* centralized assembled matrix input */
1187:     if (!mumps->myid) {
1188:       mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
1189:       if (mumps->id.ICNTL(6)>1) {
1190:         mumps->id.a = (MumpsScalar*)mumps->val;
1191:       }
1192:     }
1193:     break;
1194:   case 3:  /* distributed assembled matrix input (size>1) */
1195:     mumps->id.nz_loc = mumps->nz;
1196:     mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
1197:     if (mumps->id.ICNTL(6)>1) {
1198:       mumps->id.a_loc = (MumpsScalar*)mumps->val;
1199:     }
1200:     /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
1201:     if (!mumps->myid) {
1202:       VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);
1203:       ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
1204:     } else {
1205:       VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
1206:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
1207:     }
1208:     MatCreateVecs(A,NULL,&b);
1209:     VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
1210:     ISDestroy(&is_iden);
1211:     VecDestroy(&b);
1212:     break;
1213:   }
1214:   PetscMUMPS_c(&mumps->id);
1215:   if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1));

1217:   F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS;
1218:   F->ops->solve                 = MatSolve_MUMPS;
1219:   F->ops->solvetranspose        = MatSolve_MUMPS;
1220:   F->ops->matsolve              = MatMatSolve_MUMPS;
1221: #if defined(PETSC_USE_COMPLEX)
1222:   F->ops->getinertia = NULL;
1223: #else
1224:   F->ops->getinertia = MatGetInertia_SBAIJMUMPS;
1225: #endif
1226:   return(0);
1227: }

1231: PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer)
1232: {
1233:   PetscErrorCode    ierr;
1234:   PetscBool         iascii;
1235:   PetscViewerFormat format;
1236:   Mat_MUMPS         *mumps=(Mat_MUMPS*)A->spptr;

1239:   /* check if matrix is mumps type */
1240:   if (A->ops->solve != MatSolve_MUMPS) return(0);

1242:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1243:   if (iascii) {
1244:     PetscViewerGetFormat(viewer,&format);
1245:     if (format == PETSC_VIEWER_ASCII_INFO) {
1246:       PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");
1247:       PetscViewerASCIIPrintf(viewer,"  SYM (matrix type):                   %d \n",mumps->id.sym);
1248:       PetscViewerASCIIPrintf(viewer,"  PAR (host participation):            %d \n",mumps->id.par);
1249:       PetscViewerASCIIPrintf(viewer,"  ICNTL(1) (output for error):         %d \n",mumps->id.ICNTL(1));
1250:       PetscViewerASCIIPrintf(viewer,"  ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));
1251:       PetscViewerASCIIPrintf(viewer,"  ICNTL(3) (output for global info):   %d \n",mumps->id.ICNTL(3));
1252:       PetscViewerASCIIPrintf(viewer,"  ICNTL(4) (level of printing):        %d \n",mumps->id.ICNTL(4));
1253:       PetscViewerASCIIPrintf(viewer,"  ICNTL(5) (input mat struct):         %d \n",mumps->id.ICNTL(5));
1254:       PetscViewerASCIIPrintf(viewer,"  ICNTL(6) (matrix prescaling):        %d \n",mumps->id.ICNTL(6));
1255:       PetscViewerASCIIPrintf(viewer,"  ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));
1256:       PetscViewerASCIIPrintf(viewer,"  ICNTL(8) (scalling strategy):        %d \n",mumps->id.ICNTL(8));
1257:       PetscViewerASCIIPrintf(viewer,"  ICNTL(10) (max num of refinements):  %d \n",mumps->id.ICNTL(10));
1258:       PetscViewerASCIIPrintf(viewer,"  ICNTL(11) (error analysis):          %d \n",mumps->id.ICNTL(11));
1259:       if (mumps->id.ICNTL(11)>0) {
1260:         PetscViewerASCIIPrintf(viewer,"    RINFOG(4) (inf norm of input mat):        %g\n",mumps->id.RINFOG(4));
1261:         PetscViewerASCIIPrintf(viewer,"    RINFOG(5) (inf norm of solution):         %g\n",mumps->id.RINFOG(5));
1262:         PetscViewerASCIIPrintf(viewer,"    RINFOG(6) (inf norm of residual):         %g\n",mumps->id.RINFOG(6));
1263:         PetscViewerASCIIPrintf(viewer,"    RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));
1264:         PetscViewerASCIIPrintf(viewer,"    RINFOG(9) (error estimate):               %g \n",mumps->id.RINFOG(9));
1265:         PetscViewerASCIIPrintf(viewer,"    RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));
1266:       }
1267:       PetscViewerASCIIPrintf(viewer,"  ICNTL(12) (efficiency control):                         %d \n",mumps->id.ICNTL(12));
1268:       PetscViewerASCIIPrintf(viewer,"  ICNTL(13) (efficiency control):                         %d \n",mumps->id.ICNTL(13));
1269:       PetscViewerASCIIPrintf(viewer,"  ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));
1270:       /* ICNTL(15-17) not used */
1271:       PetscViewerASCIIPrintf(viewer,"  ICNTL(18) (input mat struct):                           %d \n",mumps->id.ICNTL(18));
1272:       PetscViewerASCIIPrintf(viewer,"  ICNTL(19) (Shur complement info):                       %d \n",mumps->id.ICNTL(19));
1273:       PetscViewerASCIIPrintf(viewer,"  ICNTL(20) (rhs sparse pattern):                         %d \n",mumps->id.ICNTL(20));
1274:       PetscViewerASCIIPrintf(viewer,"  ICNTL(21) (solution struct):                            %d \n",mumps->id.ICNTL(21));
1275:       PetscViewerASCIIPrintf(viewer,"  ICNTL(22) (in-core/out-of-core facility):               %d \n",mumps->id.ICNTL(22));
1276:       PetscViewerASCIIPrintf(viewer,"  ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));

1278:       PetscViewerASCIIPrintf(viewer,"  ICNTL(24) (detection of null pivot rows):               %d \n",mumps->id.ICNTL(24));
1279:       PetscViewerASCIIPrintf(viewer,"  ICNTL(25) (computation of a null space basis):          %d \n",mumps->id.ICNTL(25));
1280:       PetscViewerASCIIPrintf(viewer,"  ICNTL(26) (Schur options for rhs or solution):          %d \n",mumps->id.ICNTL(26));
1281:       PetscViewerASCIIPrintf(viewer,"  ICNTL(27) (experimental parameter):                     %d \n",mumps->id.ICNTL(27));
1282:       PetscViewerASCIIPrintf(viewer,"  ICNTL(28) (use parallel or sequential ordering):        %d \n",mumps->id.ICNTL(28));
1283:       PetscViewerASCIIPrintf(viewer,"  ICNTL(29) (parallel ordering):                          %d \n",mumps->id.ICNTL(29));

1285:       PetscViewerASCIIPrintf(viewer,"  ICNTL(30) (user-specified set of entries in inv(A)):    %d \n",mumps->id.ICNTL(30));
1286:       PetscViewerASCIIPrintf(viewer,"  ICNTL(31) (factors is discarded in the solve phase):    %d \n",mumps->id.ICNTL(31));
1287:       PetscViewerASCIIPrintf(viewer,"  ICNTL(33) (compute determinant):                        %d \n",mumps->id.ICNTL(33));

1289:       PetscViewerASCIIPrintf(viewer,"  CNTL(1) (relative pivoting threshold):      %g \n",mumps->id.CNTL(1));
1290:       PetscViewerASCIIPrintf(viewer,"  CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));
1291:       PetscViewerASCIIPrintf(viewer,"  CNTL(3) (absolute pivoting threshold):      %g \n",mumps->id.CNTL(3));
1292:       PetscViewerASCIIPrintf(viewer,"  CNTL(4) (value of static pivoting):         %g \n",mumps->id.CNTL(4));
1293:       PetscViewerASCIIPrintf(viewer,"  CNTL(5) (fixation for null pivots):         %g \n",mumps->id.CNTL(5));

1295:       /* infomation local to each processor */
1296:       PetscViewerASCIIPrintf(viewer, "  RINFO(1) (local estimated flops for the elimination after analysis): \n");
1297:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
1298:       PetscViewerASCIISynchronizedPrintf(viewer,"    [%d] %g \n",mumps->myid,mumps->id.RINFO(1));
1299:       PetscViewerFlush(viewer);
1300:       PetscViewerASCIIPrintf(viewer, "  RINFO(2) (local estimated flops for the assembly after factorization): \n");
1301:       PetscViewerASCIISynchronizedPrintf(viewer,"    [%d]  %g \n",mumps->myid,mumps->id.RINFO(2));
1302:       PetscViewerFlush(viewer);
1303:       PetscViewerASCIIPrintf(viewer, "  RINFO(3) (local estimated flops for the elimination after factorization): \n");
1304:       PetscViewerASCIISynchronizedPrintf(viewer,"    [%d]  %g \n",mumps->myid,mumps->id.RINFO(3));
1305:       PetscViewerFlush(viewer);

1307:       PetscViewerASCIIPrintf(viewer, "  INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");
1308:       PetscViewerASCIISynchronizedPrintf(viewer,"  [%d] %d \n",mumps->myid,mumps->id.INFO(15));
1309:       PetscViewerFlush(viewer);

1311:       PetscViewerASCIIPrintf(viewer, "  INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");
1312:       PetscViewerASCIISynchronizedPrintf(viewer,"    [%d] %d \n",mumps->myid,mumps->id.INFO(16));
1313:       PetscViewerFlush(viewer);

1315:       PetscViewerASCIIPrintf(viewer, "  INFO(23) (num of pivots eliminated on this processor after factorization): \n");
1316:       PetscViewerASCIISynchronizedPrintf(viewer,"    [%d] %d \n",mumps->myid,mumps->id.INFO(23));
1317:       PetscViewerFlush(viewer);

1319:       if (mumps->ninfo && mumps->ninfo <= 40){
1320:         PetscInt i;
1321:         for (i=0; i<mumps->ninfo; i++){
1322:           PetscViewerASCIIPrintf(viewer, "  INFO(%d): \n",mumps->info[i]);
1323:           PetscViewerASCIISynchronizedPrintf(viewer,"    [%d] %d \n",mumps->myid,mumps->id.INFO(mumps->info[i]));
1324:           PetscViewerFlush(viewer);
1325:         }
1326:       }


1329:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);

1331:       if (!mumps->myid) { /* information from the host */
1332:         PetscViewerASCIIPrintf(viewer,"  RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));
1333:         PetscViewerASCIIPrintf(viewer,"  RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));
1334:         PetscViewerASCIIPrintf(viewer,"  RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));
1335:         PetscViewerASCIIPrintf(viewer,"  (RINFOG(12) RINFOG(13))*2^INFOG(34) (determinant): (%g,%g)*(2^%d)\n",mumps->id.RINFOG(12),mumps->id.RINFOG(13),mumps->id.INFOG(34));

1337:         PetscViewerASCIIPrintf(viewer,"  INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));
1338:         PetscViewerASCIIPrintf(viewer,"  INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));
1339:         PetscViewerASCIIPrintf(viewer,"  INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));
1340:         PetscViewerASCIIPrintf(viewer,"  INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));
1341:         PetscViewerASCIIPrintf(viewer,"  INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));
1342:         PetscViewerASCIIPrintf(viewer,"  INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));
1343:         PetscViewerASCIIPrintf(viewer,"  INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));
1344:         PetscViewerASCIIPrintf(viewer,"  INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));
1345:         PetscViewerASCIIPrintf(viewer,"  INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));
1346:         PetscViewerASCIIPrintf(viewer,"  INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));
1347:         PetscViewerASCIIPrintf(viewer,"  INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));
1348:         PetscViewerASCIIPrintf(viewer,"  INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));
1349:         PetscViewerASCIIPrintf(viewer,"  INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));
1350:         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",mumps->id.INFOG(16));
1351:         PetscViewerASCIIPrintf(viewer,"  INFOG(17) (estimated size of all MUMPS internal data for factorization after analysis: sum over all processors): %d \n",mumps->id.INFOG(17));
1352:         PetscViewerASCIIPrintf(viewer,"  INFOG(18) (size of all MUMPS internal data allocated during factorization: value on the most memory consuming processor): %d \n",mumps->id.INFOG(18));
1353:         PetscViewerASCIIPrintf(viewer,"  INFOG(19) (size of all MUMPS internal data allocated during factorization: sum over all processors): %d \n",mumps->id.INFOG(19));
1354:         PetscViewerASCIIPrintf(viewer,"  INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));
1355:         PetscViewerASCIIPrintf(viewer,"  INFOG(21) (size in MB of memory effectively used during factorization - value on the most memory consuming processor): %d \n",mumps->id.INFOG(21));
1356:         PetscViewerASCIIPrintf(viewer,"  INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): %d \n",mumps->id.INFOG(22));
1357:         PetscViewerASCIIPrintf(viewer,"  INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));
1358:         PetscViewerASCIIPrintf(viewer,"  INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));
1359:         PetscViewerASCIIPrintf(viewer,"  INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));
1360:         PetscViewerASCIIPrintf(viewer,"  INFOG(28) (after factorization: number of null pivots encountered): %d\n",mumps->id.INFOG(28));
1361:         PetscViewerASCIIPrintf(viewer,"  INFOG(29) (after factorization: effective number of entries in the factors (sum over all processors)): %d\n",mumps->id.INFOG(29));
1362:         PetscViewerASCIIPrintf(viewer,"  INFOG(30, 31) (after solution: size in Mbytes of memory used during solution phase): %d, %d\n",mumps->id.INFOG(30),mumps->id.INFOG(31));
1363:         PetscViewerASCIIPrintf(viewer,"  INFOG(32) (after analysis: type of analysis done): %d\n",mumps->id.INFOG(32));
1364:         PetscViewerASCIIPrintf(viewer,"  INFOG(33) (value used for ICNTL(8)): %d\n",mumps->id.INFOG(33));
1365:         PetscViewerASCIIPrintf(viewer,"  INFOG(34) (exponent of the determinant if determinant is requested): %d\n",mumps->id.INFOG(34));
1366:       }
1367:     }
1368:   }
1369:   return(0);
1370: }

1374: PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info)
1375: {
1376:   Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr;

1379:   info->block_size        = 1.0;
1380:   info->nz_allocated      = mumps->id.INFOG(20);
1381:   info->nz_used           = mumps->id.INFOG(20);
1382:   info->nz_unneeded       = 0.0;
1383:   info->assemblies        = 0.0;
1384:   info->mallocs           = 0.0;
1385:   info->memory            = 0.0;
1386:   info->fill_ratio_given  = 0;
1387:   info->fill_ratio_needed = 0;
1388:   info->factor_mallocs    = 0;
1389:   return(0);
1390: }

1392: /* -------------------------------------------------------------------------------------------*/
1395: PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival)
1396: {
1397:   Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;

1400:   mumps->id.ICNTL(icntl) = ival;
1401:   return(0);
1402: }

1406: PetscErrorCode MatMumpsGetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt *ival)
1407: {
1408:   Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;

1411:   *ival = mumps->id.ICNTL(icntl);
1412:   return(0);
1413: }

1417: /*@
1418:   MatMumpsSetIcntl - Set MUMPS parameter ICNTL()

1420:    Logically Collective on Mat

1422:    Input Parameters:
1423: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1424: .  icntl - index of MUMPS parameter array ICNTL()
1425: -  ival - value of MUMPS ICNTL(icntl)

1427:   Options Database:
1428: .   -mat_mumps_icntl_<icntl> <ival>

1430:    Level: beginner

1432:    References: MUMPS Users' Guide

1434: .seealso: MatGetFactor()
1435: @*/
1436: PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival)
1437: {

1443:   PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
1444:   return(0);
1445: }

1449: /*@
1450:   MatMumpsGetIcntl - Get MUMPS parameter ICNTL()

1452:    Logically Collective on Mat

1454:    Input Parameters:
1455: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1456: -  icntl - index of MUMPS parameter array ICNTL()

1458:   Output Parameter:
1459: .  ival - value of MUMPS ICNTL(icntl)

1461:    Level: beginner

1463:    References: MUMPS Users' Guide

1465: .seealso: MatGetFactor()
1466: @*/
1467: PetscErrorCode MatMumpsGetIcntl(Mat F,PetscInt icntl,PetscInt *ival)
1468: {

1474:   PetscTryMethod(F,"MatMumpsGetIcntl_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));
1475:   return(0);
1476: }

1478: /* -------------------------------------------------------------------------------------------*/
1481: PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val)
1482: {
1483:   Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;

1486:   mumps->id.CNTL(icntl) = val;
1487:   return(0);
1488: }

1492: PetscErrorCode MatMumpsGetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal *val)
1493: {
1494:   Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;

1497:   *val = mumps->id.CNTL(icntl);
1498:   return(0);
1499: }

1503: /*@
1504:   MatMumpsSetCntl - Set MUMPS parameter CNTL()

1506:    Logically Collective on Mat

1508:    Input Parameters:
1509: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1510: .  icntl - index of MUMPS parameter array CNTL()
1511: -  val - value of MUMPS CNTL(icntl)

1513:   Options Database:
1514: .   -mat_mumps_cntl_<icntl> <val>

1516:    Level: beginner

1518:    References: MUMPS Users' Guide

1520: .seealso: MatGetFactor()
1521: @*/
1522: PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val)
1523: {

1529:   PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));
1530:   return(0);
1531: }

1535: /*@
1536:   MatMumpsGetCntl - Get MUMPS parameter CNTL()

1538:    Logically Collective on Mat

1540:    Input Parameters:
1541: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1542: -  icntl - index of MUMPS parameter array CNTL()

1544:   Output Parameter:
1545: .  val - value of MUMPS CNTL(icntl)

1547:    Level: beginner

1549:    References: MUMPS Users' Guide

1551: .seealso: MatGetFactor()
1552: @*/
1553: PetscErrorCode MatMumpsGetCntl(Mat F,PetscInt icntl,PetscReal *val)
1554: {

1560:   PetscTryMethod(F,"MatMumpsGetCntl_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));
1561:   return(0);
1562: }

1566: PetscErrorCode MatMumpsGetInfo_MUMPS(Mat F,PetscInt icntl,PetscInt *info)
1567: {
1568:   Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;

1571:   *info = mumps->id.INFO(icntl);
1572:   return(0);
1573: }

1577: PetscErrorCode MatMumpsGetInfog_MUMPS(Mat F,PetscInt icntl,PetscInt *infog)
1578: {
1579:   Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;

1582:   *infog = mumps->id.INFOG(icntl);
1583:   return(0);
1584: }

1588: PetscErrorCode MatMumpsGetRinfo_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfo)
1589: {
1590:   Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;

1593:   *rinfo = mumps->id.RINFO(icntl);
1594:   return(0);
1595: }

1599: PetscErrorCode MatMumpsGetRinfog_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfog)
1600: {
1601:   Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;

1604:   *rinfog = mumps->id.RINFOG(icntl);
1605:   return(0);
1606: }

1610: /*@
1611:   MatMumpsGetInfo - Get MUMPS parameter INFO()

1613:    Logically Collective on Mat

1615:    Input Parameters:
1616: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1617: -  icntl - index of MUMPS parameter array INFO()

1619:   Output Parameter:
1620: .  ival - value of MUMPS INFO(icntl)

1622:    Level: beginner

1624:    References: MUMPS Users' Guide

1626: .seealso: MatGetFactor()
1627: @*/
1628: PetscErrorCode MatMumpsGetInfo(Mat F,PetscInt icntl,PetscInt *ival)
1629: {

1634:   PetscTryMethod(F,"MatMumpsGetInfo_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));
1635:   return(0);
1636: }

1640: /*@
1641:   MatMumpsGetInfog - Get MUMPS parameter INFOG()

1643:    Logically Collective on Mat

1645:    Input Parameters:
1646: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1647: -  icntl - index of MUMPS parameter array INFOG()

1649:   Output Parameter:
1650: .  ival - value of MUMPS INFOG(icntl)

1652:    Level: beginner

1654:    References: MUMPS Users' Guide

1656: .seealso: MatGetFactor()
1657: @*/
1658: PetscErrorCode MatMumpsGetInfog(Mat F,PetscInt icntl,PetscInt *ival)
1659: {

1664:   PetscTryMethod(F,"MatMumpsGetInfog_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));
1665:   return(0);
1666: }

1670: /*@
1671:   MatMumpsGetRinfo - Get MUMPS parameter RINFO()

1673:    Logically Collective on Mat

1675:    Input Parameters:
1676: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1677: -  icntl - index of MUMPS parameter array RINFO()

1679:   Output Parameter:
1680: .  val - value of MUMPS RINFO(icntl)

1682:    Level: beginner

1684:    References: MUMPS Users' Guide

1686: .seealso: MatGetFactor()
1687: @*/
1688: PetscErrorCode MatMumpsGetRinfo(Mat F,PetscInt icntl,PetscReal *val)
1689: {

1694:   PetscTryMethod(F,"MatMumpsGetRinfo_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));
1695:   return(0);
1696: }

1700: /*@
1701:   MatMumpsGetRinfog - Get MUMPS parameter RINFOG()

1703:    Logically Collective on Mat

1705:    Input Parameters:
1706: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1707: -  icntl - index of MUMPS parameter array RINFOG()

1709:   Output Parameter:
1710: .  val - value of MUMPS RINFOG(icntl)

1712:    Level: beginner

1714:    References: MUMPS Users' Guide

1716: .seealso: MatGetFactor()
1717: @*/
1718: PetscErrorCode MatMumpsGetRinfog(Mat F,PetscInt icntl,PetscReal *val)
1719: {

1724:   PetscTryMethod(F,"MatMumpsGetRinfog_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));
1725:   return(0);
1726: }

1728: /*MC
1729:   MATSOLVERMUMPS -  A matrix type providing direct solvers (LU and Cholesky) for
1730:   distributed and sequential matrices via the external package MUMPS.

1732:   Works with MATAIJ and MATSBAIJ matrices

1734:   Options Database Keys:
1735: +  -mat_mumps_icntl_1 <6>: ICNTL(1): output stream for error messages (None)
1736: .  -mat_mumps_icntl_2 <0>: ICNTL(2): output stream for diagnostic printing, statistics, and warning (None)
1737: .  -mat_mumps_icntl_3 <0>: ICNTL(3): output stream for global information, collected on the host (None)
1738: .  -mat_mumps_icntl_4 <0>: ICNTL(4): level of printing (0 to 4) (None)
1739: .  -mat_mumps_icntl_6 <7>: ICNTL(6): permutes to a zero-free diagonal and/or scale the matrix (0 to 7) (None)
1740: .  -mat_mumps_icntl_7 <7>: ICNTL(7): computes a symmetric permutation in sequential analysis (0 to 7). 3=Scotch, 4=PORD, 5=Metis (None)
1741: .  -mat_mumps_icntl_8 <77>: ICNTL(8): scaling strategy (-2 to 8 or 77) (None)
1742: .  -mat_mumps_icntl_10 <0>: ICNTL(10): max num of refinements (None)
1743: .  -mat_mumps_icntl_11 <0>: ICNTL(11): statistics related to an error analysis (via -ksp_view) (None)
1744: .  -mat_mumps_icntl_12 <1>: ICNTL(12): an ordering strategy for symmetric matrices (0 to 3) (None)
1745: .  -mat_mumps_icntl_13 <0>: ICNTL(13): parallelism of the root node (enable ScaLAPACK) and its splitting (None)
1746: .  -mat_mumps_icntl_14 <20>: ICNTL(14): percentage increase in the estimated working space (None)
1747: .  -mat_mumps_icntl_19 <0>: ICNTL(19): computes the Schur complement (None)
1748: .  -mat_mumps_icntl_22 <0>: ICNTL(22): in-core/out-of-core factorization and solve (0 or 1) (None)
1749: .  -mat_mumps_icntl_23 <0>: ICNTL(23): max size of the working memory (MB) that can allocate per processor (None)
1750: .  -mat_mumps_icntl_24 <0>: ICNTL(24): detection of null pivot rows (0 or 1) (None)
1751: .  -mat_mumps_icntl_25 <0>: ICNTL(25): compute a solution of a deficient matrix and a null space basis (None)
1752: .  -mat_mumps_icntl_26 <0>: ICNTL(26): drives the solution phase if a Schur complement matrix (None)
1753: .  -mat_mumps_icntl_28 <1>: ICNTL(28): use 1 for sequential analysis and ictnl(7) ordering, or 2 for parallel analysis and ictnl(29) ordering (None)
1754: .  -mat_mumps_icntl_29 <0>: ICNTL(29): parallel ordering 1 = ptscotch, 2 = parmetis (None)
1755: .  -mat_mumps_icntl_30 <0>: ICNTL(30): compute user-specified set of entries in inv(A) (None)
1756: .  -mat_mumps_icntl_31 <0>: ICNTL(31): indicates which factors may be discarded during factorization (None)
1757: .  -mat_mumps_icntl_33 <0>: ICNTL(33): compute determinant (None)
1758: .  -mat_mumps_cntl_1 <0.01>: CNTL(1): relative pivoting threshold (None)
1759: .  -mat_mumps_cntl_2 <1.49012e-08>: CNTL(2): stopping criterion of refinement (None)
1760: .  -mat_mumps_cntl_3 <0>: CNTL(3): absolute pivoting threshold (None)
1761: .  -mat_mumps_cntl_4 <-1>: CNTL(4): value for static pivoting (None)
1762: -  -mat_mumps_cntl_5 <0>: CNTL(5): fixation for null pivots (None)

1764:   Level: beginner

1766: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage

1768: M*/

1772: static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type)
1773: {
1775:   *type = MATSOLVERMUMPS;
1776:   return(0);
1777: }

1779: /* MatGetFactor for Seq and MPI AIJ matrices */
1782: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F)
1783: {
1784:   Mat            B;
1786:   Mat_MUMPS      *mumps;
1787:   PetscBool      isSeqAIJ;

1790:   /* Create the factorization matrix */
1791:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
1792:   MatCreate(PetscObjectComm((PetscObject)A),&B);
1793:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1794:   MatSetType(B,((PetscObject)A)->type_name);
1795:   if (isSeqAIJ) {
1796:     MatSeqAIJSetPreallocation(B,0,NULL);
1797:   } else {
1798:     MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1799:   }

1801:   PetscNewLog(B,&mumps);

1803:   B->ops->view        = MatView_MUMPS;
1804:   B->ops->getinfo     = MatGetInfo_MUMPS;
1805:   B->ops->getdiagonal = MatGetDiagonal_MUMPS;

1807:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1808:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);
1809:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);
1810:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);
1811:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);

1813:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);
1814:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);
1815:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);
1816:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);
1817:   if (ftype == MAT_FACTOR_LU) {
1818:     B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS;
1819:     B->factortype            = MAT_FACTOR_LU;
1820:     if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij;
1821:     else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij;
1822:     mumps->sym = 0;
1823:   } else {
1824:     B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS;
1825:     B->factortype                  = MAT_FACTOR_CHOLESKY;
1826:     if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij;
1827:     else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij;
1828:     if (A->spd_set && A->spd) mumps->sym = 1;
1829:     else                      mumps->sym = 2;
1830:   }

1832:   mumps->isAIJ    = PETSC_TRUE;
1833:   mumps->Destroy  = B->ops->destroy;
1834:   B->ops->destroy = MatDestroy_MUMPS;
1835:   B->spptr        = (void*)mumps;

1837:   PetscInitializeMUMPS(A,mumps);

1839:   *F = B;
1840:   return(0);
1841: }

1843: /* MatGetFactor for Seq and MPI SBAIJ matrices */
1846: PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F)
1847: {
1848:   Mat            B;
1850:   Mat_MUMPS      *mumps;
1851:   PetscBool      isSeqSBAIJ;

1854:   if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix");
1855:   if (A->rmap->bs > 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with block size > 1 with MUMPS Cholesky, use AIJ matrix instead");
1856:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
1857:   /* Create the factorization matrix */
1858:   MatCreate(PetscObjectComm((PetscObject)A),&B);
1859:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1860:   MatSetType(B,((PetscObject)A)->type_name);
1861:   PetscNewLog(B,&mumps);
1862:   if (isSeqSBAIJ) {
1863:     MatSeqSBAIJSetPreallocation(B,1,0,NULL);

1865:     mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij;
1866:   } else {
1867:     MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);

1869:     mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij;
1870:   }

1872:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS;
1873:   B->ops->view                   = MatView_MUMPS;
1874:   B->ops->getdiagonal            = MatGetDiagonal_MUMPS;

1876:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1877:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);
1878:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);
1879:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);
1880:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);

1882:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);
1883:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);
1884:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);
1885:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);

1887:   B->factortype = MAT_FACTOR_CHOLESKY;
1888:   if (A->spd_set && A->spd) mumps->sym = 1;
1889:   else                      mumps->sym = 2;

1891:   mumps->isAIJ    = PETSC_FALSE;
1892:   mumps->Destroy  = B->ops->destroy;
1893:   B->ops->destroy = MatDestroy_MUMPS;
1894:   B->spptr        = (void*)mumps;

1896:   PetscInitializeMUMPS(A,mumps);

1898:   *F = B;
1899:   return(0);
1900: }

1904: PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F)
1905: {
1906:   Mat            B;
1908:   Mat_MUMPS      *mumps;
1909:   PetscBool      isSeqBAIJ;

1912:   /* Create the factorization matrix */
1913:   PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);
1914:   MatCreate(PetscObjectComm((PetscObject)A),&B);
1915:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1916:   MatSetType(B,((PetscObject)A)->type_name);
1917:   if (isSeqBAIJ) {
1918:     MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);
1919:   } else {
1920:     MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);
1921:   }

1923:   PetscNewLog(B,&mumps);
1924:   if (ftype == MAT_FACTOR_LU) {
1925:     B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS;
1926:     B->factortype            = MAT_FACTOR_LU;
1927:     if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij;
1928:     else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij;
1929:     mumps->sym = 0;
1930:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n");

1932:   B->ops->view        = MatView_MUMPS;
1933:   B->ops->getdiagonal = MatGetDiagonal_MUMPS;

1935:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1936:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);
1937:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);
1938:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);
1939:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);

1941:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);
1942:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);
1943:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);
1944:   PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);

1946:   mumps->isAIJ    = PETSC_TRUE;
1947:   mumps->Destroy  = B->ops->destroy;
1948:   B->ops->destroy = MatDestroy_MUMPS;
1949:   B->spptr        = (void*)mumps;

1951:   PetscInitializeMUMPS(A,mumps);

1953:   *F = B;
1954:   return(0);
1955: }

1957: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*);
1958: PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat,MatFactorType,Mat*);
1959: PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*);

1963: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MUMPS(void)
1964: {

1968:   MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ,         MAT_FACTOR_LU,MatGetFactor_aij_mumps);
1969:   MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ,         MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);
1970:   MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ,         MAT_FACTOR_LU,MatGetFactor_baij_mumps);
1971:   MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ,         MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);
1972:   MatSolverPackageRegister(MATSOLVERMUMPS,MATMPISBAIJ,         MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);
1973:   MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ,         MAT_FACTOR_LU,MatGetFactor_aij_mumps);
1974:   MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ,         MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);
1975:   MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ,         MAT_FACTOR_LU,MatGetFactor_baij_mumps);
1976:   MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ,         MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);
1977:   MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQSBAIJ,         MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);
1978:   return(0);
1979: }