Actual source code: mumps.c
petsc-3.4.2 2013-07-02
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));
669: if (nneg) {
670: if (!mumps->myid) {
671: *nneg = mumps->id.INFOG(12);
672: }
673: MPI_Bcast(nneg,1,MPI_INT,0,mumps->comm_mumps);
674: }
675: if (nzero) *nzero = 0;
676: if (npos) *npos = F->rmap->N - (*nneg);
677: return(0);
678: }
679: #endif /* !defined(PETSC_USE_COMPLEX) */
683: PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info)
684: {
685: Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr;
687: Mat F_diag;
688: PetscBool isMPIAIJ;
691: (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);
693: /* numerical factorization phase */
694: /*-------------------------------*/
695: mumps->id.job = JOB_FACTNUMERIC;
696: if (!mumps->id.ICNTL(18)) {
697: if (!mumps->myid) {
698: #if defined(PETSC_USE_COMPLEX)
699: #if defined(PETSC_USE_REAL_SINGLE)
700: mumps->id.a = (mumps_complex*)mumps->val;
701: #else
702: mumps->id.a = (mumps_double_complex*)mumps->val;
703: #endif
704: #else
705: mumps->id.a = mumps->val;
706: #endif
707: }
708: } else {
709: #if defined(PETSC_USE_COMPLEX)
710: #if defined(PETSC_USE_REAL_SINGLE)
711: mumps->id.a_loc = (mumps_complex*)mumps->val;
712: #else
713: mumps->id.a_loc = (mumps_double_complex*)mumps->val;
714: #endif
715: #else
716: mumps->id.a_loc = mumps->val;
717: #endif
718: }
719: PetscMUMPS_c(&mumps->id);
720: if (mumps->id.INFOG(1) < 0) {
721: if (mumps->id.INFO(1) == -13) 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));
722: 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));
723: }
724: 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));
726: if (mumps->size > 1) {
727: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
728: if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A;
729: else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A;
730: F_diag->assembled = PETSC_TRUE;
731: if (mumps->scat_sol) {
732: VecScatterDestroy(&mumps->scat_sol);
733: PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);
734: VecDestroy(&mumps->x_seq);
735: }
736: }
737: (F)->assembled = PETSC_TRUE;
738: mumps->matstruc = SAME_NONZERO_PATTERN;
739: mumps->CleanUpMUMPS = PETSC_TRUE;
741: if (mumps->size > 1) {
742: /* distributed solution */
743: if (!mumps->scat_sol) {
744: /* Create x_seq=sol_loc for repeated use */
745: PetscInt lsol_loc;
746: PetscScalar *sol_loc;
748: lsol_loc = mumps->id.INFO(23); /* length of sol_loc */
750: PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&mumps->id.isol_loc);
752: mumps->id.lsol_loc = lsol_loc;
753: #if defined(PETSC_USE_COMPLEX)
754: #if defined(PETSC_USE_REAL_SINGLE)
755: mumps->id.sol_loc = (mumps_complex*)sol_loc;
756: #else
757: mumps->id.sol_loc = (mumps_double_complex*)sol_loc;
758: #endif
759: #else
760: mumps->id.sol_loc = sol_loc;
761: #endif
762: VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);
763: }
764: }
765: return(0);
766: }
768: /* Sets MUMPS options from the options database */
771: PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A)
772: {
773: Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr;
775: PetscInt icntl;
776: PetscBool flg;
779: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");
780: PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);
781: if (flg) mumps->id.ICNTL(1) = icntl;
782: PetscOptionsInt("-mat_mumps_icntl_2","ICNTL(2): output stream for diagnostic printing, statistics, and warning","None",mumps->id.ICNTL(2),&icntl,&flg);
783: if (flg) mumps->id.ICNTL(2) = icntl;
784: PetscOptionsInt("-mat_mumps_icntl_3","ICNTL(3): output stream for global information, collected on the host","None",mumps->id.ICNTL(3),&icntl,&flg);
785: if (flg) mumps->id.ICNTL(3) = icntl;
787: PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);
788: if (flg) mumps->id.ICNTL(4) = icntl;
789: if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */
791: PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): permuting and/or scaling the matrix (0 to 7)","None",mumps->id.ICNTL(6),&icntl,&flg);
792: if (flg) mumps->id.ICNTL(6) = icntl;
794: 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);
795: if (flg) {
796: 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");
797: else mumps->id.ICNTL(7) = icntl;
798: }
800: PetscOptionsInt("-mat_mumps_icntl_8","ICNTL(8): scaling strategy (-2 to 8 or 77)","None",mumps->id.ICNTL(8),&mumps->id.ICNTL(8),NULL);
801: PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);
802: 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);
803: 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);
804: PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control: with or without ScaLAPACK","None",mumps->id.ICNTL(13),&mumps->id.ICNTL(13),NULL);
805: PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): percentage of estimated workspace increase","None",mumps->id.ICNTL(14),&mumps->id.ICNTL(14),NULL);
806: PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);
808: 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);
809: 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);
810: 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);
811: if (mumps->id.ICNTL(24)) {
812: mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */
813: }
815: PetscOptionsInt("-mat_mumps_icntl_25","ICNTL(25): computation of a null space basis","None",mumps->id.ICNTL(25),&mumps->id.ICNTL(25),NULL);
816: 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);
817: PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);
818: 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);
819: PetscOptionsInt("-mat_mumps_icntl_29","ICNTL(29): parallel ordering 1 = ptscotch 2 = parmetis","None",mumps->id.ICNTL(29),&mumps->id.ICNTL(29),NULL);
820: 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);
821: 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);
822: PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);
824: PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);
825: PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);
826: PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);
827: PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);
828: PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);
830: PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL);
831: PetscOptionsEnd();
832: return(0);
833: }
837: PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps)
838: {
842: MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid);
843: MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);
844: MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));
846: mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps);
848: mumps->id.job = JOB_INIT;
849: mumps->id.par = 1; /* host participates factorizaton and solve */
850: mumps->id.sym = mumps->sym;
851: PetscMUMPS_c(&mumps->id);
853: mumps->CleanUpMUMPS = PETSC_FALSE;
854: mumps->scat_rhs = NULL;
855: mumps->scat_sol = NULL;
857: /* set PETSc-MUMPS default options - override MUMPS default */
858: mumps->id.ICNTL(3) = 0;
859: mumps->id.ICNTL(4) = 0;
860: if (mumps->size == 1) {
861: mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */
862: } else {
863: mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */
864: mumps->id.ICNTL(21) = 1; /* distributed solution */
865: }
866: return(0);
867: }
869: /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */
872: PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
873: {
874: Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr;
876: Vec b;
877: IS is_iden;
878: const PetscInt M = A->rmap->N;
881: mumps->matstruc = DIFFERENT_NONZERO_PATTERN;
883: /* Set MUMPS options from the options database */
884: PetscSetMUMPSFromOptions(F,A);
886: (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);
888: /* analysis phase */
889: /*----------------*/
890: mumps->id.job = JOB_FACTSYMBOLIC;
891: mumps->id.n = M;
892: switch (mumps->id.ICNTL(18)) {
893: case 0: /* centralized assembled matrix input */
894: if (!mumps->myid) {
895: mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
896: if (mumps->id.ICNTL(6)>1) {
897: #if defined(PETSC_USE_COMPLEX)
898: #if defined(PETSC_USE_REAL_SINGLE)
899: mumps->id.a = (mumps_complex*)mumps->val;
900: #else
901: mumps->id.a = (mumps_double_complex*)mumps->val;
902: #endif
903: #else
904: mumps->id.a = mumps->val;
905: #endif
906: }
907: if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */
908: /*
909: PetscBool flag;
910: ISEqual(r,c,&flag);
911: if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm");
912: ISView(r,PETSC_VIEWER_STDOUT_SELF);
913: */
914: if (!mumps->myid) {
915: const PetscInt *idx;
916: PetscInt i,*perm_in;
918: PetscMalloc(M*sizeof(PetscInt),&perm_in);
919: ISGetIndices(r,&idx);
921: mumps->id.perm_in = perm_in;
922: for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */
923: ISRestoreIndices(r,&idx);
924: }
925: }
926: }
927: break;
928: case 3: /* distributed assembled matrix input (size>1) */
929: mumps->id.nz_loc = mumps->nz;
930: mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
931: if (mumps->id.ICNTL(6)>1) {
932: #if defined(PETSC_USE_COMPLEX)
933: #if defined(PETSC_USE_REAL_SINGLE)
934: mumps->id.a_loc = (mumps_complex*)mumps->val;
935: #else
936: mumps->id.a_loc = (mumps_double_complex*)mumps->val;
937: #endif
938: #else
939: mumps->id.a_loc = mumps->val;
940: #endif
941: }
942: /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
943: if (!mumps->myid) {
944: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);
945: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
946: } else {
947: VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
948: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
949: }
950: VecCreate(PetscObjectComm((PetscObject)A),&b);
951: VecSetSizes(b,A->rmap->n,PETSC_DECIDE);
952: VecSetFromOptions(b);
954: VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
955: ISDestroy(&is_iden);
956: VecDestroy(&b);
957: break;
958: }
959: PetscMUMPS_c(&mumps->id);
960: 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));
962: F->ops->lufactornumeric = MatFactorNumeric_MUMPS;
963: F->ops->solve = MatSolve_MUMPS;
964: F->ops->solvetranspose = MatSolveTranspose_MUMPS;
965: F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */
966: return(0);
967: }
969: /* Note the Petsc r and c permutations are ignored */
972: PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
973: {
974: Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr;
976: Vec b;
977: IS is_iden;
978: const PetscInt M = A->rmap->N;
981: mumps->matstruc = DIFFERENT_NONZERO_PATTERN;
983: /* Set MUMPS options from the options database */
984: PetscSetMUMPSFromOptions(F,A);
986: (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);
988: /* analysis phase */
989: /*----------------*/
990: mumps->id.job = JOB_FACTSYMBOLIC;
991: mumps->id.n = M;
992: switch (mumps->id.ICNTL(18)) {
993: case 0: /* centralized assembled matrix input */
994: if (!mumps->myid) {
995: mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
996: if (mumps->id.ICNTL(6)>1) {
997: #if defined(PETSC_USE_COMPLEX)
998: #if defined(PETSC_USE_REAL_SINGLE)
999: mumps->id.a = (mumps_complex*)mumps->val;
1000: #else
1001: mumps->id.a = (mumps_double_complex*)mumps->val;
1002: #endif
1003: #else
1004: mumps->id.a = mumps->val;
1005: #endif
1006: }
1007: }
1008: break;
1009: case 3: /* distributed assembled matrix input (size>1) */
1010: mumps->id.nz_loc = mumps->nz;
1011: mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
1012: if (mumps->id.ICNTL(6)>1) {
1013: #if defined(PETSC_USE_COMPLEX)
1014: #if defined(PETSC_USE_REAL_SINGLE)
1015: mumps->id.a_loc = (mumps_complex*)mumps->val;
1016: #else
1017: mumps->id.a_loc = (mumps_double_complex*)mumps->val;
1018: #endif
1019: #else
1020: mumps->id.a_loc = mumps->val;
1021: #endif
1022: }
1023: /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
1024: if (!mumps->myid) {
1025: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);
1026: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
1027: } else {
1028: VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
1029: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
1030: }
1031: VecCreate(PetscObjectComm((PetscObject)A),&b);
1032: VecSetSizes(b,A->rmap->n,PETSC_DECIDE);
1033: VecSetFromOptions(b);
1035: VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
1036: ISDestroy(&is_iden);
1037: VecDestroy(&b);
1038: break;
1039: }
1040: PetscMUMPS_c(&mumps->id);
1041: 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));
1043: F->ops->lufactornumeric = MatFactorNumeric_MUMPS;
1044: F->ops->solve = MatSolve_MUMPS;
1045: F->ops->solvetranspose = MatSolveTranspose_MUMPS;
1046: return(0);
1047: }
1049: /* Note the Petsc r permutation and factor info are ignored */
1052: PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info)
1053: {
1054: Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr;
1056: Vec b;
1057: IS is_iden;
1058: const PetscInt M = A->rmap->N;
1061: mumps->matstruc = DIFFERENT_NONZERO_PATTERN;
1063: /* Set MUMPS options from the options database */
1064: PetscSetMUMPSFromOptions(F,A);
1066: (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);
1068: /* analysis phase */
1069: /*----------------*/
1070: mumps->id.job = JOB_FACTSYMBOLIC;
1071: mumps->id.n = M;
1072: switch (mumps->id.ICNTL(18)) {
1073: case 0: /* centralized assembled matrix input */
1074: if (!mumps->myid) {
1075: mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn;
1076: if (mumps->id.ICNTL(6)>1) {
1077: #if defined(PETSC_USE_COMPLEX)
1078: #if defined(PETSC_USE_REAL_SINGLE)
1079: mumps->id.a = (mumps_complex*)mumps->val;
1080: #else
1081: mumps->id.a = (mumps_double_complex*)mumps->val;
1082: #endif
1083: #else
1084: mumps->id.a = mumps->val;
1085: #endif
1086: }
1087: }
1088: break;
1089: case 3: /* distributed assembled matrix input (size>1) */
1090: mumps->id.nz_loc = mumps->nz;
1091: mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn;
1092: if (mumps->id.ICNTL(6)>1) {
1093: #if defined(PETSC_USE_COMPLEX)
1094: #if defined(PETSC_USE_REAL_SINGLE)
1095: mumps->id.a_loc = (mumps_complex*)mumps->val;
1096: #else
1097: mumps->id.a_loc = (mumps_double_complex*)mumps->val;
1098: #endif
1099: #else
1100: mumps->id.a_loc = mumps->val;
1101: #endif
1102: }
1103: /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
1104: if (!mumps->myid) {
1105: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);
1106: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
1107: } else {
1108: VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);
1109: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);
1110: }
1111: VecCreate(PetscObjectComm((PetscObject)A),&b);
1112: VecSetSizes(b,A->rmap->n,PETSC_DECIDE);
1113: VecSetFromOptions(b);
1115: VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);
1116: ISDestroy(&is_iden);
1117: VecDestroy(&b);
1118: break;
1119: }
1120: PetscMUMPS_c(&mumps->id);
1121: 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));
1123: F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS;
1124: F->ops->solve = MatSolve_MUMPS;
1125: F->ops->solvetranspose = MatSolve_MUMPS;
1126: F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */
1127: #if !defined(PETSC_USE_COMPLEX)
1128: F->ops->getinertia = MatGetInertia_SBAIJMUMPS;
1129: #else
1130: F->ops->getinertia = NULL;
1131: #endif
1132: return(0);
1133: }
1137: PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer)
1138: {
1139: PetscErrorCode ierr;
1140: PetscBool iascii;
1141: PetscViewerFormat format;
1142: Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr;
1145: /* check if matrix is mumps type */
1146: if (A->ops->solve != MatSolve_MUMPS) return(0);
1148: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1149: if (iascii) {
1150: PetscViewerGetFormat(viewer,&format);
1151: if (format == PETSC_VIEWER_ASCII_INFO) {
1152: PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");
1153: PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);
1154: PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);
1155: PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));
1156: PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));
1157: PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));
1158: PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));
1159: PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));
1160: PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));
1161: PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));
1162: PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));
1163: PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));
1164: PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));
1165: if (mumps->id.ICNTL(11)>0) {
1166: PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));
1167: PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));
1168: PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));
1169: PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));
1170: PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));
1171: PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));
1172: }
1173: PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));
1174: PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));
1175: PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));
1176: /* ICNTL(15-17) not used */
1177: PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));
1178: PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));
1179: PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));
1180: PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));
1181: PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));
1182: PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));
1184: PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));
1185: PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));
1186: PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));
1187: PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));
1188: PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));
1189: PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));
1191: PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));
1192: PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));
1193: PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));
1195: PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));
1196: PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));
1197: PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));
1198: PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));
1199: PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));
1201: /* infomation local to each processor */
1202: PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");
1203: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
1204: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));
1205: PetscViewerFlush(viewer);
1206: PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");
1207: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));
1208: PetscViewerFlush(viewer);
1209: PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");
1210: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));
1211: PetscViewerFlush(viewer);
1213: PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");
1214: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));
1215: PetscViewerFlush(viewer);
1217: PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");
1218: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));
1219: PetscViewerFlush(viewer);
1221: PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");
1222: PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));
1223: PetscViewerFlush(viewer);
1224: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
1226: if (!mumps->myid) { /* information from the host */
1227: PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));
1228: PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));
1229: PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));
1230: 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));
1232: PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));
1233: PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));
1234: PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));
1235: PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));
1236: PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));
1237: PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));
1238: PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));
1239: PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));
1240: PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));
1241: PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));
1242: PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));
1243: PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));
1244: PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));
1245: 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));
1246: 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));
1247: 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));
1248: PetscViewerASCIIPrintf(viewer," INFOG(19) (size of all MUMPS internal data allocated during factorization: sum over all processors): %d \n",mumps->id.INFOG(19));
1249: PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));
1250: 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));
1251: PetscViewerASCIIPrintf(viewer," INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): %d \n",mumps->id.INFOG(22));
1252: PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));
1253: PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));
1254: PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));
1255: }
1256: }
1257: }
1258: return(0);
1259: }
1263: PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info)
1264: {
1265: Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr;
1268: info->block_size = 1.0;
1269: info->nz_allocated = mumps->id.INFOG(20);
1270: info->nz_used = mumps->id.INFOG(20);
1271: info->nz_unneeded = 0.0;
1272: info->assemblies = 0.0;
1273: info->mallocs = 0.0;
1274: info->memory = 0.0;
1275: info->fill_ratio_given = 0;
1276: info->fill_ratio_needed = 0;
1277: info->factor_mallocs = 0;
1278: return(0);
1279: }
1281: /* -------------------------------------------------------------------------------------------*/
1284: PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival)
1285: {
1286: Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;
1289: mumps->id.ICNTL(icntl) = ival;
1290: return(0);
1291: }
1295: /*@
1296: MatMumpsSetIcntl - Set MUMPS parameter ICNTL()
1298: Logically Collective on Mat
1300: Input Parameters:
1301: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1302: . icntl - index of MUMPS parameter array ICNTL()
1303: - ival - value of MUMPS ICNTL(icntl)
1305: Options Database:
1306: . -mat_mumps_icntl_<icntl> <ival>
1308: Level: beginner
1310: References: MUMPS Users' Guide
1312: .seealso: MatGetFactor()
1313: @*/
1314: PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival)
1315: {
1321: PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
1322: return(0);
1323: }
1325: /* -------------------------------------------------------------------------------------------*/
1328: PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val)
1329: {
1330: Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr;
1333: mumps->id.CNTL(icntl) = val;
1334: return(0);
1335: }
1339: /*@
1340: MatMumpsSetCntl - Set MUMPS parameter CNTL()
1342: Logically Collective on Mat
1344: Input Parameters:
1345: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface
1346: . icntl - index of MUMPS parameter array CNTL()
1347: - val - value of MUMPS CNTL(icntl)
1349: Options Database:
1350: . -mat_mumps_cntl_<icntl> <val>
1352: Level: beginner
1354: References: MUMPS Users' Guide
1356: .seealso: MatGetFactor()
1357: @*/
1358: PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val)
1359: {
1365: PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));
1366: return(0);
1367: }
1369: /*MC
1370: MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for
1371: distributed and sequential matrices via the external package MUMPS.
1373: Works with MATAIJ and MATSBAIJ matrices
1375: Options Database Keys:
1376: + -mat_mumps_icntl_4 <0,...,4> - print level
1377: . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide)
1378: . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec)
1379: . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T
1380: . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements
1381: . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view
1382: . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide)
1383: . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide)
1384: . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide)
1385: . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide)
1386: . -mat_mumps_cntl_1 <delta> - relative pivoting threshold
1387: . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement
1388: - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold
1390: Level: beginner
1392: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
1394: M*/
1398: static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type)
1399: {
1401: *type = MATSOLVERMUMPS;
1402: return(0);
1403: }
1405: /* MatGetFactor for Seq and MPI AIJ matrices */
1408: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F)
1409: {
1410: Mat B;
1412: Mat_MUMPS *mumps;
1413: PetscBool isSeqAIJ;
1416: /* Create the factorization matrix */
1417: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
1418: MatCreate(PetscObjectComm((PetscObject)A),&B);
1419: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1420: MatSetType(B,((PetscObject)A)->type_name);
1421: if (isSeqAIJ) {
1422: MatSeqAIJSetPreallocation(B,0,NULL);
1423: } else {
1424: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
1425: }
1427: PetscNewLog(B,Mat_MUMPS,&mumps);
1429: B->ops->view = MatView_MUMPS;
1430: B->ops->getinfo = MatGetInfo_MUMPS;
1432: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1433: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);
1434: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);
1435: if (ftype == MAT_FACTOR_LU) {
1436: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS;
1437: B->factortype = MAT_FACTOR_LU;
1438: if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij;
1439: else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij;
1440: mumps->sym = 0;
1441: } else {
1442: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS;
1443: B->factortype = MAT_FACTOR_CHOLESKY;
1444: if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij;
1445: else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij;
1446: if (A->spd_set && A->spd) mumps->sym = 1;
1447: else mumps->sym = 2;
1448: }
1450: mumps->isAIJ = PETSC_TRUE;
1451: mumps->Destroy = B->ops->destroy;
1452: B->ops->destroy = MatDestroy_MUMPS;
1453: B->spptr = (void*)mumps;
1455: PetscInitializeMUMPS(A,mumps);
1457: *F = B;
1458: return(0);
1459: }
1461: /* MatGetFactor for Seq and MPI SBAIJ matrices */
1464: PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F)
1465: {
1466: Mat B;
1468: Mat_MUMPS *mumps;
1469: PetscBool isSeqSBAIJ;
1472: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix");
1473: 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");
1474: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
1475: /* Create the factorization matrix */
1476: MatCreate(PetscObjectComm((PetscObject)A),&B);
1477: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1478: MatSetType(B,((PetscObject)A)->type_name);
1479: PetscNewLog(B,Mat_MUMPS,&mumps);
1480: if (isSeqSBAIJ) {
1481: MatSeqSBAIJSetPreallocation(B,1,0,NULL);
1483: mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij;
1484: } else {
1485: MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);
1487: mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij;
1488: }
1490: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS;
1491: B->ops->view = MatView_MUMPS;
1493: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1494: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl);
1495: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl);
1497: B->factortype = MAT_FACTOR_CHOLESKY;
1498: if (A->spd_set && A->spd) mumps->sym = 1;
1499: else mumps->sym = 2;
1501: mumps->isAIJ = PETSC_FALSE;
1502: mumps->Destroy = B->ops->destroy;
1503: B->ops->destroy = MatDestroy_MUMPS;
1504: B->spptr = (void*)mumps;
1506: PetscInitializeMUMPS(A,mumps);
1508: *F = B;
1509: return(0);
1510: }
1514: PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F)
1515: {
1516: Mat B;
1518: Mat_MUMPS *mumps;
1519: PetscBool isSeqBAIJ;
1522: /* Create the factorization matrix */
1523: PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);
1524: MatCreate(PetscObjectComm((PetscObject)A),&B);
1525: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1526: MatSetType(B,((PetscObject)A)->type_name);
1527: if (isSeqBAIJ) {
1528: MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);
1529: } else {
1530: MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);
1531: }
1533: PetscNewLog(B,Mat_MUMPS,&mumps);
1534: if (ftype == MAT_FACTOR_LU) {
1535: B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS;
1536: B->factortype = MAT_FACTOR_LU;
1537: if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij;
1538: else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij;
1539: mumps->sym = 0;
1540: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n");
1542: B->ops->view = MatView_MUMPS;
1544: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);
1545: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);
1546: PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);
1548: mumps->isAIJ = PETSC_TRUE;
1549: mumps->Destroy = B->ops->destroy;
1550: B->ops->destroy = MatDestroy_MUMPS;
1551: B->spptr = (void*)mumps;
1553: PetscInitializeMUMPS(A,mumps);
1555: *F = B;
1556: return(0);
1557: }