Actual source code: mumps.c
petsc-3.4.5 2014-06-29
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
46: /* macros s.t. indices match MUMPS documentation */
47: #define ICNTL(I) icntl[(I)-1]
48: #define CNTL(I) cntl[(I)-1]
49: #define INFOG(I) infog[(I)-1]
50: #define INFO(I) info[(I)-1]
51: #define RINFOG(I) rinfog[(I)-1]
52: #define RINFO(I) rinfo[(I)-1]
54: typedef struct {
55: #if defined(PETSC_USE_COMPLEX)
56: #if defined(PETSC_USE_REAL_SINGLE)
57: CMUMPS_STRUC_C id;
58: #else
59: ZMUMPS_STRUC_C id;
60: #endif
61: #else
62: #if defined(PETSC_USE_REAL_SINGLE)
63: SMUMPS_STRUC_C id;
64: #else
65: DMUMPS_STRUC_C id;
66: #endif
67: #endif
69: MatStructure matstruc;
70: PetscMPIInt myid,size;
71: PetscInt *irn,*jcn,nz,sym;
72: PetscScalar *val;
73: MPI_Comm comm_mumps;
74: VecScatter scat_rhs, scat_sol;
75: PetscBool isAIJ,CleanUpMUMPS;
76: Vec b_seq,x_seq;
77: PetscInt ICNTL9_pre; /* check if ICNTL(9) is changed from previous MatSolve */
79: PetscErrorCode (*Destroy)(Mat);
80: PetscErrorCode (*ConvertToTriples)(Mat, int, MatReuse, int*, int**, int**, PetscScalar**);
81: } Mat_MUMPS;
83: extern PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*);
86: /* MatConvertToTriples_A_B */
87: /*convert Petsc matrix to triples: row[nz], col[nz], val[nz] */
88: /*
89: input:
90: A - matrix in aij,baij or sbaij (bs=1) format
91: shift - 0: C style output triple; 1: Fortran style output triple.
92: reuse - MAT_INITIAL_MATRIX: spaces are allocated and values are set for the triple
93: MAT_REUSE_MATRIX: only the values in v array are updated
94: output:
95: nnz - dim of r, c, and v (number of local nonzero entries of A)
96: r, c, v - row and col index, matrix values (matrix triples)
97: */
101: PetscErrorCode MatConvertToTriples_seqaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
102: {
103: const PetscInt *ai,*aj,*ajj,M=A->rmap->n;
104: PetscInt nz,rnz,i,j;
106: PetscInt *row,*col;
107: Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data;
110: *v=aa->a;
111: if (reuse == MAT_INITIAL_MATRIX) {
112: nz = aa->nz;
113: ai = aa->i;
114: aj = aa->j;
115: *nnz = nz;
116: PetscMalloc(2*nz*sizeof(PetscInt), &row);
117: col = row + nz;
119: nz = 0;
120: for (i=0; i<M; i++) {
121: rnz = ai[i+1] - ai[i];
122: ajj = aj + ai[i];
123: for (j=0; j<rnz; j++) {
124: row[nz] = i+shift; col[nz++] = ajj[j] + shift;
125: }
126: }
127: *r = row; *c = col;
128: }
129: return(0);
130: }
134: PetscErrorCode MatConvertToTriples_seqbaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
135: {
136: Mat_SeqBAIJ *aa=(Mat_SeqBAIJ*)A->data;
137: const PetscInt *ai,*aj,*ajj,bs=A->rmap->bs,bs2=aa->bs2,M=A->rmap->N/bs;
138: PetscInt nz,idx=0,rnz,i,j,k,m;
140: PetscInt *row,*col;
143: *v = aa->a;
144: if (reuse == MAT_INITIAL_MATRIX) {
145: ai = aa->i; aj = aa->j;
146: nz = bs2*aa->nz;
147: *nnz = nz;
148: PetscMalloc(2*nz*sizeof(PetscInt), &row);
149: col = row + nz;
151: for (i=0; i<M; i++) {
152: ajj = aj + ai[i];
153: rnz = ai[i+1] - ai[i];
154: for (k=0; k<rnz; k++) {
155: for (j=0; j<bs; j++) {
156: for (m=0; m<bs; m++) {
157: row[idx] = i*bs + m + shift;
158: col[idx++] = bs*(ajj[k]) + j + shift;
159: }
160: }
161: }
162: }
163: *r = row; *c = col;
164: }
165: return(0);
166: }
170: PetscErrorCode MatConvertToTriples_seqsbaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
171: {
172: const PetscInt *ai, *aj,*ajj,M=A->rmap->n;
173: PetscInt nz,rnz,i,j;
175: PetscInt *row,*col;
176: Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data;
179: *v = aa->a;
180: if (reuse == MAT_INITIAL_MATRIX) {
181: nz = aa->nz;
182: ai = aa->i;
183: aj = aa->j;
184: *v = aa->a;
185: *nnz = nz;
186: PetscMalloc(2*nz*sizeof(PetscInt), &row);
187: col = row + nz;
189: nz = 0;
190: for (i=0; i<M; i++) {
191: rnz = ai[i+1] - ai[i];
192: ajj = aj + ai[i];
193: for (j=0; j<rnz; j++) {
194: row[nz] = i+shift; col[nz++] = ajj[j] + shift;
195: }
196: }
197: *r = row; *c = col;
198: }
199: return(0);
200: }
204: PetscErrorCode MatConvertToTriples_seqaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
205: {
206: const PetscInt *ai,*aj,*ajj,*adiag,M=A->rmap->n;
207: PetscInt nz,rnz,i,j;
208: const PetscScalar *av,*v1;
209: PetscScalar *val;
210: PetscErrorCode ierr;
211: PetscInt *row,*col;
212: Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data;
215: ai =aa->i; aj=aa->j;av=aa->a;
216: adiag=aa->diag;
217: if (reuse == MAT_INITIAL_MATRIX) {
218: nz = M + (aa->nz-M)/2;
219: *nnz = nz;
220: PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
221: col = row + nz;
222: val = (PetscScalar*)(col + nz);
224: nz = 0;
225: for (i=0; i<M; i++) {
226: rnz = ai[i+1] - adiag[i];
227: ajj = aj + adiag[i];
228: v1 = av + adiag[i];
229: for (j=0; j<rnz; j++) {
230: row[nz] = i+shift; col[nz] = ajj[j] + shift; val[nz++] = v1[j];
231: }
232: }
233: *r = row; *c = col; *v = val;
234: } else {
235: nz = 0; val = *v;
236: for (i=0; i <M; i++) {
237: rnz = ai[i+1] - adiag[i];
238: ajj = aj + adiag[i];
239: v1 = av + adiag[i];
240: for (j=0; j<rnz; j++) {
241: val[nz++] = v1[j];
242: }
243: }
244: }
245: return(0);
246: }
250: PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
251: {
252: const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
253: PetscErrorCode ierr;
254: PetscInt rstart,nz,i,j,jj,irow,countA,countB;
255: PetscInt *row,*col;
256: const PetscScalar *av, *bv,*v1,*v2;
257: PetscScalar *val;
258: Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data;
259: Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)(mat->A)->data;
260: Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data;
263: ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart;
264: av=aa->a; bv=bb->a;
266: garray = mat->garray;
268: if (reuse == MAT_INITIAL_MATRIX) {
269: nz = aa->nz + bb->nz;
270: *nnz = nz;
271: PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
272: col = row + nz;
273: val = (PetscScalar*)(col + nz);
275: *r = row; *c = col; *v = val;
276: } else {
277: row = *r; col = *c; val = *v;
278: }
280: jj = 0; irow = rstart;
281: for (i=0; i<m; i++) {
282: ajj = aj + ai[i]; /* ptr to the beginning of this row */
283: countA = ai[i+1] - ai[i];
284: countB = bi[i+1] - bi[i];
285: bjj = bj + bi[i];
286: v1 = av + ai[i];
287: v2 = bv + bi[i];
289: /* A-part */
290: for (j=0; j<countA; j++) {
291: if (reuse == MAT_INITIAL_MATRIX) {
292: row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift;
293: }
294: val[jj++] = v1[j];
295: }
297: /* B-part */
298: for (j=0; j < countB; j++) {
299: if (reuse == MAT_INITIAL_MATRIX) {
300: row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift;
301: }
302: val[jj++] = v2[j];
303: }
304: irow++;
305: }
306: return(0);
307: }
311: PetscErrorCode MatConvertToTriples_mpiaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
312: {
313: const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
314: PetscErrorCode ierr;
315: PetscInt rstart,nz,i,j,jj,irow,countA,countB;
316: PetscInt *row,*col;
317: const PetscScalar *av, *bv,*v1,*v2;
318: PetscScalar *val;
319: Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data;
320: Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(mat->A)->data;
321: Mat_SeqAIJ *bb = (Mat_SeqAIJ*)(mat->B)->data;
324: ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart;
325: av=aa->a; bv=bb->a;
327: garray = mat->garray;
329: if (reuse == MAT_INITIAL_MATRIX) {
330: nz = aa->nz + bb->nz;
331: *nnz = nz;
332: PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
333: col = row + nz;
334: val = (PetscScalar*)(col + nz);
336: *r = row; *c = col; *v = val;
337: } else {
338: row = *r; col = *c; val = *v;
339: }
341: jj = 0; irow = rstart;
342: for (i=0; i<m; i++) {
343: ajj = aj + ai[i]; /* ptr to the beginning of this row */
344: countA = ai[i+1] - ai[i];
345: countB = bi[i+1] - bi[i];
346: bjj = bj + bi[i];
347: v1 = av + ai[i];
348: v2 = bv + bi[i];
350: /* A-part */
351: for (j=0; j<countA; j++) {
352: if (reuse == MAT_INITIAL_MATRIX) {
353: row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift;
354: }
355: val[jj++] = v1[j];
356: }
358: /* B-part */
359: for (j=0; j < countB; j++) {
360: if (reuse == MAT_INITIAL_MATRIX) {
361: row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift;
362: }
363: val[jj++] = v2[j];
364: }
365: irow++;
366: }
367: return(0);
368: }
372: PetscErrorCode MatConvertToTriples_mpibaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
373: {
374: Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data;
375: Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)(mat->A)->data;
376: Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data;
377: const PetscInt *ai = aa->i, *bi = bb->i, *aj = aa->j, *bj = bb->j,*ajj, *bjj;
378: const PetscInt *garray = mat->garray,mbs=mat->mbs,rstart=A->rmap->rstart;
379: const PetscInt bs = A->rmap->bs,bs2=mat->bs2;
380: PetscErrorCode ierr;
381: PetscInt nz,i,j,k,n,jj,irow,countA,countB,idx;
382: PetscInt *row,*col;
383: const PetscScalar *av=aa->a, *bv=bb->a,*v1,*v2;
384: PetscScalar *val;
387: if (reuse == MAT_INITIAL_MATRIX) {
388: nz = bs2*(aa->nz + bb->nz);
389: *nnz = nz;
390: PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
391: col = row + nz;
392: val = (PetscScalar*)(col + nz);
394: *r = row; *c = col; *v = val;
395: } else {
396: row = *r; col = *c; val = *v;
397: }
399: jj = 0; irow = rstart;
400: for (i=0; i<mbs; i++) {
401: countA = ai[i+1] - ai[i];
402: countB = bi[i+1] - bi[i];
403: ajj = aj + ai[i];
404: bjj = bj + bi[i];
405: v1 = av + bs2*ai[i];
406: v2 = bv + bs2*bi[i];
408: idx = 0;
409: /* A-part */
410: for (k=0; k<countA; k++) {
411: for (j=0; j<bs; j++) {
412: for (n=0; n<bs; n++) {
413: if (reuse == MAT_INITIAL_MATRIX) {
414: row[jj] = irow + n + shift;
415: col[jj] = rstart + bs*ajj[k] + j + shift;
416: }
417: val[jj++] = v1[idx++];
418: }
419: }
420: }
422: idx = 0;
423: /* B-part */
424: for (k=0; k<countB; k++) {
425: for (j=0; j<bs; j++) {
426: for (n=0; n<bs; n++) {
427: if (reuse == MAT_INITIAL_MATRIX) {
428: row[jj] = irow + n + shift;
429: col[jj] = bs*garray[bjj[k]] + j + shift;
430: }
431: val[jj++] = v2[idx++];
432: }
433: }
434: }
435: irow += bs;
436: }
437: return(0);
438: }
442: PetscErrorCode MatConvertToTriples_mpiaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v)
443: {
444: const PetscInt *ai, *aj,*adiag, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
445: PetscErrorCode ierr;
446: PetscInt rstart,nz,nza,nzb,i,j,jj,irow,countA,countB;
447: PetscInt *row,*col;
448: const PetscScalar *av, *bv,*v1,*v2;
449: PetscScalar *val;
450: Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data;
451: Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(mat->A)->data;
452: Mat_SeqAIJ *bb =(Mat_SeqAIJ*)(mat->B)->data;
455: ai=aa->i; aj=aa->j; adiag=aa->diag;
456: bi=bb->i; bj=bb->j; garray = mat->garray;
457: av=aa->a; bv=bb->a;
459: rstart = A->rmap->rstart;
461: if (reuse == MAT_INITIAL_MATRIX) {
462: nza = 0; /* num of upper triangular entries in mat->A, including diagonals */
463: nzb = 0; /* num of upper triangular entries in mat->B */
464: for (i=0; i<m; i++) {
465: nza += (ai[i+1] - adiag[i]);
466: countB = bi[i+1] - bi[i];
467: bjj = bj + bi[i];
468: for (j=0; j<countB; j++) {
469: if (garray[bjj[j]] > rstart) nzb++;
470: }
471: }
473: nz = nza + nzb; /* total nz of upper triangular part of mat */
474: *nnz = nz;
475: PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);
476: col = row + nz;
477: val = (PetscScalar*)(col + nz);
479: *r = row; *c = col; *v = val;
480: } else {
481: row = *r; col = *c; val = *v;
482: }
484: jj = 0; irow = rstart;
485: for (i=0; i<m; i++) {
486: ajj = aj + adiag[i]; /* ptr to the beginning of the diagonal of this row */
487: v1 = av + adiag[i];
488: countA = ai[i+1] - adiag[i];
489: countB = bi[i+1] - bi[i];
490: bjj = bj + bi[i];
491: v2 = bv + bi[i];
493: /* A-part */
494: for (j=0; j<countA; j++) {
495: if (reuse == MAT_INITIAL_MATRIX) {
496: row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift;
497: }
498: val[jj++] = v1[j];
499: }
501: /* B-part */
502: for (j=0; j < countB; j++) {
503: if (garray[bjj[j]] > rstart) {
504: if (reuse == MAT_INITIAL_MATRIX) {
505: row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift;
506: }
507: val[jj++] = v2[j];
508: }
509: }
510: irow++;
511: }
512: return(0);
513: }
517: PetscErrorCode MatDestroy_MUMPS(Mat A)
518: {
519: Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr;
523: if (mumps->CleanUpMUMPS) {
524: /* Terminate instance, deallocate memories */
525: PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);
526: VecScatterDestroy(&mumps->scat_rhs);
527: VecDestroy(&mumps->b_seq);
528: VecScatterDestroy(&mumps->scat_sol);
529: VecDestroy(&mumps->x_seq);
530: PetscFree(mumps->id.perm_in);
531: PetscFree(mumps->irn);
533: mumps->id.job = JOB_END;
534: PetscMUMPS_c(&mumps->id);
535: MPI_Comm_free(&(mumps->comm_mumps));
536: }
537: if (mumps->Destroy) {
538: (mumps->Destroy)(A);
539: }
540: PetscFree(A->spptr);
542: /* clear composed functions */
543: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
544: PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetIcntl_C",NULL);
545: PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetCntl_C",NULL);
546: return(0);
547: }
551: PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x)
552: {
553: Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr;
554: PetscScalar *array;
555: Vec b_seq;
556: IS is_iden,is_petsc;
558: PetscInt i;
561: mumps->id.nrhs = 1;
562: b_seq = mumps->b_seq;
563: if (mumps->size > 1) {
564: /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */
565: VecScatterBegin(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);
566: VecScatterEnd(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);
567: if (!mumps->myid) {VecGetArray(b_seq,&array);}
568: } else { /* size == 1 */
569: VecCopy(b,x);
570: VecGetArray(x,&array);
571: }
572: if (!mumps->myid) { /* define rhs on the host */
573: mumps->id.nrhs = 1;
574: #if defined(PETSC_USE_COMPLEX)
575: #if defined(PETSC_USE_REAL_SINGLE)
576: mumps->id.rhs = (mumps_complex*)array;
577: #else
578: mumps->id.rhs = (mumps_double_complex*)array;
579: #endif
580: #else
581: mumps->id.rhs = array;
582: #endif
583: }
585: /* solve phase */
586: /*-------------*/
587: mumps->id.job = JOB_SOLVE;
588: PetscMUMPS_c(&mumps->id);
589: 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));
591: if (mumps->size > 1) { /* convert mumps distributed solution to petsc mpi x */
592: if (mumps->scat_sol && mumps->ICNTL9_pre != mumps->id.ICNTL(9)) {
593: /* when id.ICNTL(9) changes, the contents of lsol_loc may change (not its size, lsol_loc), recreates scat_sol */
594: VecScatterDestroy(&mumps->scat_sol);
595: }
596: if (!mumps->scat_sol) { /* create scatter scat_sol */
597: ISCreateStride(PETSC_COMM_SELF,mumps->id.lsol_loc,0,1,&is_iden); /* from */
598: for (i=0; i<mumps->id.lsol_loc; i++) {
599: mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */
600: }
601: ISCreateGeneral(PETSC_COMM_SELF,mumps->id.lsol_loc,mumps->id.isol_loc,PETSC_COPY_VALUES,&is_petsc); /* to */
602: VecScatterCreate(mumps->x_seq,is_iden,x,is_petsc,&mumps->scat_sol);
603: ISDestroy(&is_iden);
604: ISDestroy(&is_petsc);
606: mumps->ICNTL9_pre = mumps->id.ICNTL(9); /* save current value of id.ICNTL(9) */
607: }
609: VecScatterBegin(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
610: VecScatterEnd(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
611: }
612: return(0);
613: }
617: PetscErrorCode MatSolveTranspose_MUMPS(Mat A,Vec b,Vec x)
618: {
619: Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr;
623: mumps->id.ICNTL(9) = 0;
625: MatSolve_MUMPS(A,b,x);
627: mumps->id.ICNTL(9) = 1;
628: return(0);
629: }
633: PetscErrorCode MatMatSolve_MUMPS(Mat A,Mat B,Mat X)
634: {
636: PetscBool flg;
639: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
640: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
641: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
642: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
643: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_MUMPS() is not implemented yet");
644: return(0);
645: }
647: #if !defined(PETSC_USE_COMPLEX)
648: /*
649: input:
650: F: numeric factor
651: output:
652: nneg: total number of negative pivots
653: nzero: 0
654: npos: (global dimension of F) - nneg
655: */
659: PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos)
660: {
661: Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;
663: PetscMPIInt size;
666: MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);
667: /* 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 */
668: 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));
670: if (nneg) *nneg = mumps->id.INFOG(12);
671: if (nzero || npos) {
672: 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");
673: if (nzero) *nzero = mumps->id.INFOG(28);
674: if (npos) *npos = F->rmap->N - (mumps->id.INFOG(12) + mumps->id.INFOG(28));
675: }
676: return(0);
677: }
678: #endif /* !defined(PETSC_USE_COMPLEX) */
682: PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info)
683: {
684: Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr;
686: Mat F_diag;
687: PetscBool isMPIAIJ;
690: (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);
692: /* numerical factorization phase */
693: /*-------------------------------*/
694: mumps->id.job = JOB_FACTNUMERIC;
695: if (!mumps->id.ICNTL(18)) {
696: if (!mumps->myid) {
697: #if defined(PETSC_USE_COMPLEX)
698: #if defined(PETSC_USE_REAL_SINGLE)
699: mumps->id.a = (mumps_complex*)mumps->val;
700: #else
701: mumps->id.a = (mumps_double_complex*)mumps->val;
702: #endif
703: #else
704: mumps->id.a = mumps->val;
705: #endif
706: }
707: } else {
708: #if defined(PETSC_USE_COMPLEX)
709: #if defined(PETSC_USE_REAL_SINGLE)
710: mumps->id.a_loc = (mumps_complex*)mumps->val;
711: #else
712: mumps->id.a_loc = (mumps_double_complex*)mumps->val;
713: #endif
714: #else
715: mumps->id.a_loc = mumps->val;
716: #endif
717: }
718: PetscMUMPS_c(&mumps->id);
719: if (mumps->id.INFOG(1) < 0) {
720: if (mumps->id.INFO(1) == -13) {
721: if (mumps->id.INFO(2) < 0) {
722: 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));
723: } else {
724: 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));
725: }
726: } 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));
727: }
728: 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));
730: if (mumps->size > 1) {
731: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
732: if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A;
733: else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A;
734: F_diag->assembled = PETSC_TRUE;
735: if (mumps->scat_sol) {
736: VecScatterDestroy(&mumps->scat_sol);
737: PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);
738: VecDestroy(&mumps->x_seq);
739: }
740: }
741: (F)->assembled = PETSC_TRUE;
742: mumps->matstruc = SAME_NONZERO_PATTERN;
743: mumps->CleanUpMUMPS = PETSC_TRUE;
745: if (mumps->size > 1) {
746: /* distributed solution */
747: if (!mumps->scat_sol) {
748: /* Create x_seq=sol_loc for repeated use */
749: PetscInt lsol_loc;
750: PetscScalar *sol_loc;
752: lsol_loc = mumps->id.INFO(23); /* length of sol_loc */
754: PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&mumps->id.isol_loc);
756: mumps->id.lsol_loc = lsol_loc;
757: #if defined(PETSC_USE_COMPLEX)
758: #if defined(PETSC_USE_REAL_SINGLE)
759: mumps->id.sol_loc = (mumps_complex*)sol_loc;
760: #else
761: mumps->id.sol_loc = (mumps_double_complex*)sol_loc;
762: #endif
763: #else
764: mumps->id.sol_loc = sol_loc;
765: #endif
766: VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);
767: }
768: }
769: return(0);
770: }
772: /* Sets MUMPS options from the options database */
775: PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A)
776: {
777: Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr;
779: PetscInt icntl;
780: PetscBool flg;
783: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");
784: PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);
785: if (flg) mumps->id.ICNTL(1) = icntl;
786: PetscOptionsInt("-mat_mumps_icntl_2","ICNTL(2): output stream for diagnostic printing, statistics, and warning","None",mumps->id.ICNTL(2),&icntl,&flg);
787: if (flg) mumps->id.ICNTL(2) = icntl;
788: PetscOptionsInt("-mat_mumps_icntl_3","ICNTL(3): output stream for global information, collected on the host","None",mumps->id.ICNTL(3),&icntl,&flg);
789: if (flg) mumps->id.ICNTL(3) = icntl;
791: PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);
792: if (flg) mumps->id.ICNTL(4) = icntl;
793: if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */
795: PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): permuting and/or scaling the matrix (0 to 7)","None",mumps->id.ICNTL(6),&icntl,&flg);
796: if (flg) mumps->id.ICNTL(6) = icntl;
798: PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7). 3=Scotch, 4=PORD, 5=Metis","None",mumps->id.ICNTL(7),&icntl,&flg);
799: if (flg) {
800: 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");
801: else mumps->id.ICNTL(7) = icntl;
802: }
804: PetscOptionsInt("-mat_mumps_icntl_8","ICNTL(8): scaling strategy (-2 to 8 or 77)","None",mumps->id.ICNTL(8),&mumps->id.ICNTL(8),NULL);
805: PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);
806: PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): statistics related to the linear system solved (via -ksp_view)","None",mumps->id.ICNTL(11),&mumps->id.ICNTL(11),NULL);
807: PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control: defines the ordering strategy with scaling constraints (0 to 3)","None",mumps->id.ICNTL(12),&mumps->id.ICNTL(12),NULL);
808: PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control: with or without ScaLAPACK","None",mumps->id.ICNTL(13),&mumps->id.ICNTL(13),NULL);
809: PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): percentage of estimated workspace increase","None",mumps->id.ICNTL(14),&mumps->id.ICNTL(14),NULL);
810: PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);
812: PetscOptionsInt("-mat_mumps_icntl_22","ICNTL(22): in-core/out-of-core facility (0 or 1)","None",mumps->id.ICNTL(22),&mumps->id.ICNTL(22),NULL);
813: 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);
814: 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);
815: if (mumps->id.ICNTL(24)) {
816: mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */
817: }
819: PetscOptionsInt("-mat_mumps_icntl_25","ICNTL(25): computation of a null space basis","None",mumps->id.ICNTL(25),&mumps->id.ICNTL(25),NULL);
820: PetscOptionsInt("-mat_mumps_icntl_26","ICNTL(26): Schur options for right-hand side or solution vector","None",mumps->id.ICNTL(26),&mumps->id.ICNTL(26),NULL);
821: PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);
822: 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);
823: PetscOptionsInt("-mat_mumps_icntl_29","ICNTL(29): parallel ordering 1 = ptscotch 2 = parmetis","None",mumps->id.ICNTL(29),&mumps->id.ICNTL(29),NULL);
824: 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);
825: PetscOptionsInt("-mat_mumps_icntl_31","ICNTL(31): factors can be discarded in the solve phase","None",mumps->id.ICNTL(31),&mumps->id.ICNTL(31),NULL);
826: PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);
828: PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);
829: PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);
830: PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);
831: PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);
832: PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);
834: PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL);
835: PetscOptionsEnd();
836: return(0);
837: }
841: PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps)
842: {
846: MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid);
847: MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);
848: MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));
850: mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps);
852: mumps->id.job = JOB_INIT;
853: mumps->id.par = 1; /* host participates factorizaton and solve */
854: mumps->id.sym = mumps->sym;
855: PetscMUMPS_c(&mumps->id);
857: mumps->CleanUpMUMPS = PETSC_FALSE;
858: mumps->scat_rhs = NULL;
859: mumps->scat_sol = NULL;
861: /* set PETSc-MUMPS default options - override MUMPS default */
862: mumps->id.ICNTL(3) = 0;
863: mumps->id.ICNTL(4) = 0;
864: if (mumps->size == 1) {
865: mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */
866: } else {
867: mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */
868: mumps->id.ICNTL(21) = 1; /* distributed solution */
869: }
870: return(0);
871: }
873: /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */
876: PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
877: {
878: Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr;
880: Vec b;
881: IS is_iden;
882: const PetscInt M = A->rmap->N;
885: mumps->matstruc = DIFFERENT_NONZERO_PATTERN;
887: /* Set MUMPS options from the options database */
888: PetscSetMUMPSFromOptions(F,A);
890: (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);
892: /* analysis phase */
893: /*----------------*/
894: mumps->id.job = JOB_FACTSYMBOLIC;
895: mumps->id.n = M;
896: switch (mumps->id.ICNTL(18)) {
897: case 0: /* centralized assembled matrix input */
898: if (!mumps->myid) {
899: mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
900: if (mumps->id.ICNTL(6)>1) {
901: #if defined(PETSC_USE_COMPLEX)
902: #if defined(PETSC_USE_REAL_SINGLE)
903: mumps->id.a = (mumps_complex*)mumps->val;
904: #else
905: mumps->id.a = (mumps_double_complex*)mumps->val;
906: #endif
907: #else
908: mumps->id.a = mumps->val;
909: #endif
910: }
911: if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */
912: /*
913: PetscBool flag;
914: ISEqual(r,c,&flag);
915: if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm");
916: ISView(r,PETSC_VIEWER_STDOUT_SELF);
917: */
918: if (!mumps->myid) {
919: const PetscInt *idx;
920: PetscInt i,*perm_in;
922: PetscMalloc(M*sizeof(PetscInt),&perm_in);
923: ISGetIndices(r,&idx);
925: mumps->id.perm_in = perm_in;
926: for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */
927: ISRestoreIndices(r,&idx);
928: }
929: }
930: }
931: break;
932: case 3: /* distributed assembled matrix input (size>1) */
933: mumps->id.nz_loc = mumps->nz;
934: mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
935: if (mumps->id.ICNTL(6)>1) {
936: #if defined(PETSC_USE_COMPLEX)
937: #if defined(PETSC_USE_REAL_SINGLE)
938: mumps->id.a_loc = (mumps_complex*)mumps->val;
939: #else
940: mumps->id.a_loc = (mumps_double_complex*)mumps->val;
941: #endif
942: #else
943: mumps->id.a_loc = mumps->val;
944: #endif
945: }
946: /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
947: if (!mumps->myid) {
948: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);
949: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
950: } else {
951: VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
952: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
953: }
954: VecCreate(PetscObjectComm((PetscObject)A),&b);
955: VecSetSizes(b,A->rmap->n,PETSC_DECIDE);
956: VecSetFromOptions(b);
958: VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
959: ISDestroy(&is_iden);
960: VecDestroy(&b);
961: break;
962: }
963: PetscMUMPS_c(&mumps->id);
964: 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));
966: F->ops->lufactornumeric = MatFactorNumeric_MUMPS;
967: F->ops->solve = MatSolve_MUMPS;
968: F->ops->solvetranspose = MatSolveTranspose_MUMPS;
969: F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */
970: return(0);
971: }
973: /* Note the Petsc r and c permutations are ignored */
976: PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
977: {
978: Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr;
980: Vec b;
981: IS is_iden;
982: const PetscInt M = A->rmap->N;
985: mumps->matstruc = DIFFERENT_NONZERO_PATTERN;
987: /* Set MUMPS options from the options database */
988: PetscSetMUMPSFromOptions(F,A);
990: (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);
992: /* analysis phase */
993: /*----------------*/
994: mumps->id.job = JOB_FACTSYMBOLIC;
995: mumps->id.n = M;
996: switch (mumps->id.ICNTL(18)) {
997: case 0: /* centralized assembled matrix input */
998: if (!mumps->myid) {
999: mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
1000: if (mumps->id.ICNTL(6)>1) {
1001: #if defined(PETSC_USE_COMPLEX)
1002: #if defined(PETSC_USE_REAL_SINGLE)
1003: mumps->id.a = (mumps_complex*)mumps->val;
1004: #else
1005: mumps->id.a = (mumps_double_complex*)mumps->val;
1006: #endif
1007: #else
1008: mumps->id.a = mumps->val;
1009: #endif
1010: }
1011: }
1012: break;
1013: case 3: /* distributed assembled matrix input (size>1) */
1014: mumps->id.nz_loc = mumps->nz;
1015: mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
1016: if (mumps->id.ICNTL(6)>1) {
1017: #if defined(PETSC_USE_COMPLEX)
1018: #if defined(PETSC_USE_REAL_SINGLE)
1019: mumps->id.a_loc = (mumps_complex*)mumps->val;
1020: #else
1021: mumps->id.a_loc = (mumps_double_complex*)mumps->val;
1022: #endif
1023: #else
1024: mumps->id.a_loc = mumps->val;
1025: #endif
1026: }
1027: /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
1028: if (!mumps->myid) {
1029: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);
1030: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
1031: } else {
1032: VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
1033: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
1034: }
1035: VecCreate(PetscObjectComm((PetscObject)A),&b);
1036: VecSetSizes(b,A->rmap->n,PETSC_DECIDE);
1037: VecSetFromOptions(b);
1039: VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
1040: ISDestroy(&is_iden);
1041: VecDestroy(&b);
1042: break;
1043: }
1044: PetscMUMPS_c(&mumps->id);
1045: 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));
1047: F->ops->lufactornumeric = MatFactorNumeric_MUMPS;
1048: F->ops->solve = MatSolve_MUMPS;
1049: F->ops->solvetranspose = MatSolveTranspose_MUMPS;
1050: return(0);
1051: }
1053: /* Note the Petsc r permutation and factor info are ignored */
1056: PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info)
1057: {
1058: Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr;
1060: Vec b;
1061: IS is_iden;
1062: const PetscInt M = A->rmap->N;
1065: mumps->matstruc = DIFFERENT_NONZERO_PATTERN;
1067: /* Set MUMPS options from the options database */
1068: PetscSetMUMPSFromOptions(F,A);
1070: (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);
1072: /* analysis phase */
1073: /*----------------*/
1074: mumps->id.job = JOB_FACTSYMBOLIC;
1075: mumps->id.n = M;
1076: switch (mumps->id.ICNTL(18)) {
1077: case 0: /* centralized assembled matrix input */
1078: if (!mumps->myid) {
1079: mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
1080: if (mumps->id.ICNTL(6)>1) {
1081: #if defined(PETSC_USE_COMPLEX)
1082: #if defined(PETSC_USE_REAL_SINGLE)
1083: mumps->id.a = (mumps_complex*)mumps->val;
1084: #else
1085: mumps->id.a = (mumps_double_complex*)mumps->val;
1086: #endif
1087: #else
1088: mumps->id.a = mumps->val;
1089: #endif
1090: }
1091: }
1092: break;
1093: case 3: /* distributed assembled matrix input (size>1) */
1094: mumps->id.nz_loc = mumps->nz;
1095: mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
1096: if (mumps->id.ICNTL(6)>1) {
1097: #if defined(PETSC_USE_COMPLEX)
1098: #if defined(PETSC_USE_REAL_SINGLE)
1099: mumps->id.a_loc = (mumps_complex*)mumps->val;
1100: #else
1101: mumps->id.a_loc = (mumps_double_complex*)mumps->val;
1102: #endif
1103: #else
1104: mumps->id.a_loc = mumps->val;
1105: #endif
1106: }
1107: /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
1108: if (!mumps->myid) {
1109: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);
1110: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
1111: } else {
1112: VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
1113: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
1114: }
1115: VecCreate(PetscObjectComm((PetscObject)A),&b);
1116: VecSetSizes(b,A->rmap->n,PETSC_DECIDE);
1117: VecSetFromOptions(b);
1119: VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
1120: ISDestroy(&is_iden);
1121: VecDestroy(&b);
1122: break;
1123: }
1124: PetscMUMPS_c(&mumps->id);
1125: 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));
1127: F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS;
1128: F->ops->solve = MatSolve_MUMPS;
1129: F->ops->solvetranspose = MatSolve_MUMPS;
1130: F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */
1131: #if !defined(PETSC_USE_COMPLEX)
1132: F->ops->getinertia = MatGetInertia_SBAIJMUMPS;
1133: #else
1134: F->ops->getinertia = NULL;
1135: #endif
1136: return(0);
1137: }
1141: PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer)
1142: {
1143: PetscErrorCode ierr;
1144: PetscBool iascii;
1145: PetscViewerFormat format;
1146: Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr;
1149: /* check if matrix is mumps type */
1150: if (A->ops->solve != MatSolve_MUMPS) return(0);
1152: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1153: if (iascii) {
1154: PetscViewerGetFormat(viewer,&format);
1155: if (format == PETSC_VIEWER_ASCII_INFO) {
1156: PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");
1157: PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);
1158: PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);
1159: PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));
1160: PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));
1161: PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));
1162: PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));
1163: PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));
1164: PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));
1165: PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));
1166: PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));
1167: PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));
1168: PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));
1169: if (mumps->id.ICNTL(11)>0) {
1170: PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));
1171: PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));
1172: PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));
1173: PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));
1174: PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));
1175: PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));
1176: }
1177: PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));
1178: PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));
1179: PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));
1180: /* ICNTL(15-17) not used */
1181: PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));
1182: PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));
1183: PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));
1184: PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));
1185: PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));
1186: PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));
1188: PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));
1189: PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));
1190: PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));
1191: PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));
1192: PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));
1193: PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));
1195: PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));
1196: PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));
1197: PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));
1199: PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));
1200: PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));
1201: PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));
1202: PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));
1203: PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));
1205: /* infomation local to each processor */
1206: PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");
1207: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
1208: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));
1209: PetscViewerFlush(viewer);
1210: PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");
1211: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));
1212: PetscViewerFlush(viewer);
1213: PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");
1214: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));
1215: PetscViewerFlush(viewer);
1217: PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");
1218: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));
1219: PetscViewerFlush(viewer);
1221: PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");
1222: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));
1223: PetscViewerFlush(viewer);
1225: PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");
1226: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));
1227: PetscViewerFlush(viewer);
1228: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
1230: if (!mumps->myid) { /* information from the host */
1231: PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));
1232: PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));
1233: PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));
1234: 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));
1236: PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));
1237: PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));
1238: PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));
1239: PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));
1240: PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));
1241: PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));
1242: PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));
1243: PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));
1244: PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));
1245: PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));
1246: PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));
1247: PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));
1248: PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));
1249: 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));
1250: 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));
1251: 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));
1252: PetscViewerASCIIPrintf(viewer," INFOG(19) (size of all MUMPS internal data allocated during factorization: sum over all processors): %d \n",mumps->id.INFOG(19));
1253: PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));
1254: 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));
1255: PetscViewerASCIIPrintf(viewer," INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): %d \n",mumps->id.INFOG(22));
1256: PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));
1257: PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));
1258: PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));
1259: }
1260: }
1261: }
1262: return(0);
1263: }
1267: PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info)
1268: {
1269: Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr;
1272: info->block_size = 1.0;
1273: info->nz_allocated = mumps->id.INFOG(20);
1274: info->nz_used = mumps->id.INFOG(20);
1275: info->nz_unneeded = 0.0;
1276: info->assemblies = 0.0;
1277: info->mallocs = 0.0;
1278: info->memory = 0.0;
1279: info->fill_ratio_given = 0;
1280: info->fill_ratio_needed = 0;
1281: info->factor_mallocs = 0;
1282: return(0);
1283: }
1285: /* -------------------------------------------------------------------------------------------*/
1288: PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival)
1289: {
1290: Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;
1293: mumps->id.ICNTL(icntl) = ival;
1294: return(0);
1295: }
1299: /*@
1300: MatMumpsSetIcntl - Set MUMPS parameter ICNTL()
1302: Logically Collective on Mat
1304: Input Parameters:
1305: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1306: . icntl - index of MUMPS parameter array ICNTL()
1307: - ival - value of MUMPS ICNTL(icntl)
1309: Options Database:
1310: . -mat_mumps_icntl_<icntl> <ival>
1312: Level: beginner
1314: References: MUMPS Users' Guide
1316: .seealso: MatGetFactor()
1317: @*/
1318: PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival)
1319: {
1325: PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
1326: return(0);
1327: }
1329: /* -------------------------------------------------------------------------------------------*/
1332: PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val)
1333: {
1334: Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;
1337: mumps->id.CNTL(icntl) = val;
1338: return(0);
1339: }
1343: /*@
1344: MatMumpsSetCntl - Set MUMPS parameter CNTL()
1346: Logically Collective on Mat
1348: Input Parameters:
1349: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1350: . icntl - index of MUMPS parameter array CNTL()
1351: - val - value of MUMPS CNTL(icntl)
1353: Options Database:
1354: . -mat_mumps_cntl_<icntl> <val>
1356: Level: beginner
1358: References: MUMPS Users' Guide
1360: .seealso: MatGetFactor()
1361: @*/
1362: PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val)
1363: {
1369: PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));
1370: return(0);
1371: }
1373: /*MC
1374: MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for
1375: distributed and sequential matrices via the external package MUMPS.
1377: Works with MATAIJ and MATSBAIJ matrices
1379: Options Database Keys:
1380: + -mat_mumps_icntl_4 <0,...,4> - print level
1381: . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide)
1382: . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec)
1383: . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T
1384: . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements
1385: . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view
1386: . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide)
1387: . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide)
1388: . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide)
1389: . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide)
1390: . -mat_mumps_cntl_1 <delta> - relative pivoting threshold
1391: . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement
1392: - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold
1394: Level: beginner
1396: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
1398: M*/
1402: static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type)
1403: {
1405: *type = MATSOLVERMUMPS;
1406: return(0);
1407: }
1409: /* MatGetFactor for Seq and MPI AIJ matrices */
1412: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F)
1413: {
1414: Mat B;
1416: Mat_MUMPS *mumps;
1417: PetscBool isSeqAIJ;
1420: /* Create the factorization matrix */
1421: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
1422: MatCreate(PetscObjectComm((PetscObject)A),&B);
1423: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1424: MatSetType(B,((PetscObject)A)->type_name);
1425: if (isSeqAIJ) {
1426: MatSeqAIJSetPreallocation(B,0,NULL);
1427: } else {
1428: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1429: }
1431: PetscNewLog(B,Mat_MUMPS,&mumps);
1433: B->ops->view = MatView_MUMPS;
1434: B->ops->getinfo = MatGetInfo_MUMPS;
1436: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1437: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);
1438: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);
1439: if (ftype == MAT_FACTOR_LU) {
1440: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS;
1441: B->factortype = MAT_FACTOR_LU;
1442: if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij;
1443: else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij;
1444: mumps->sym = 0;
1445: } else {
1446: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS;
1447: B->factortype = MAT_FACTOR_CHOLESKY;
1448: if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij;
1449: else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij;
1450: if (A->spd_set && A->spd) mumps->sym = 1;
1451: else mumps->sym = 2;
1452: }
1454: mumps->isAIJ = PETSC_TRUE;
1455: mumps->Destroy = B->ops->destroy;
1456: B->ops->destroy = MatDestroy_MUMPS;
1457: B->spptr = (void*)mumps;
1459: PetscInitializeMUMPS(A,mumps);
1461: *F = B;
1462: return(0);
1463: }
1465: /* MatGetFactor for Seq and MPI SBAIJ matrices */
1468: PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F)
1469: {
1470: Mat B;
1472: Mat_MUMPS *mumps;
1473: PetscBool isSeqSBAIJ;
1476: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix");
1477: 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");
1478: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
1479: /* Create the factorization matrix */
1480: MatCreate(PetscObjectComm((PetscObject)A),&B);
1481: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1482: MatSetType(B,((PetscObject)A)->type_name);
1483: PetscNewLog(B,Mat_MUMPS,&mumps);
1484: if (isSeqSBAIJ) {
1485: MatSeqSBAIJSetPreallocation(B,1,0,NULL);
1487: mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij;
1488: } else {
1489: MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);
1491: mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij;
1492: }
1494: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS;
1495: B->ops->view = MatView_MUMPS;
1497: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1498: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl);
1499: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl);
1501: B->factortype = MAT_FACTOR_CHOLESKY;
1502: if (A->spd_set && A->spd) mumps->sym = 1;
1503: else mumps->sym = 2;
1505: mumps->isAIJ = PETSC_FALSE;
1506: mumps->Destroy = B->ops->destroy;
1507: B->ops->destroy = MatDestroy_MUMPS;
1508: B->spptr = (void*)mumps;
1510: PetscInitializeMUMPS(A,mumps);
1512: *F = B;
1513: return(0);
1514: }
1518: PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F)
1519: {
1520: Mat B;
1522: Mat_MUMPS *mumps;
1523: PetscBool isSeqBAIJ;
1526: /* Create the factorization matrix */
1527: PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);
1528: MatCreate(PetscObjectComm((PetscObject)A),&B);
1529: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1530: MatSetType(B,((PetscObject)A)->type_name);
1531: if (isSeqBAIJ) {
1532: MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);
1533: } else {
1534: MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);
1535: }
1537: PetscNewLog(B,Mat_MUMPS,&mumps);
1538: if (ftype == MAT_FACTOR_LU) {
1539: B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS;
1540: B->factortype = MAT_FACTOR_LU;
1541: if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij;
1542: else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij;
1543: mumps->sym = 0;
1544: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n");
1546: B->ops->view = MatView_MUMPS;
1548: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1549: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);
1550: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);
1552: mumps->isAIJ = PETSC_TRUE;
1553: mumps->Destroy = B->ops->destroy;
1554: B->ops->destroy = MatDestroy_MUMPS;
1555: B->spptr = (void*)mumps;
1557: PetscInitializeMUMPS(A,mumps);
1559: *F = B;
1560: return(0);
1561: }