Actual source code: pastix.c
petsc-3.4.5 2014-06-29
1: /*
2: Provides an interface to the PaStiX sparse solver
3: */
4: #include <../src/mat/impls/aij/seq/aij.h>
5: #include <../src/mat/impls/aij/mpi/mpiaij.h>
6: #include <../src/mat/impls/sbaij/seq/sbaij.h>
7: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
9: #if defined(PETSC_USE_COMPLEX) && defined(__cplusplus)
10: #define _H_COMPLEX
11: #endif
13: EXTERN_C_BEGIN
14: #include <pastix.h>
15: EXTERN_C_END
17: #if defined(PETSC_USE_COMPLEX)
18: #if defined(PETSC_USE_REAL_SINGLE)
19: #define PASTIX_CALL c_pastix
20: #define PASTIX_CHECKMATRIX c_pastix_checkMatrix
21: #define PastixScalar COMPLEX
22: #else
23: #define PASTIX_CALL z_pastix
24: #define PASTIX_CHECKMATRIX z_pastix_checkMatrix
25: #define PastixScalar DCOMPLEX
26: #endif
28: #else /* PETSC_USE_COMPLEX */
30: #if defined(PETSC_USE_REAL_SINGLE)
31: #define PASTIX_CALL s_pastix
32: #define PASTIX_CHECKMATRIX s_pastix_checkMatrix
33: #define PastixScalar float
34: #else
35: #define PASTIX_CALL d_pastix
36: #define PASTIX_CHECKMATRIX d_pastix_checkMatrix
37: #define PastixScalar double
38: #endif
40: #endif /* PETSC_USE_COMPLEX */
42: typedef struct Mat_Pastix_ {
43: pastix_data_t *pastix_data; /* Pastix data storage structure */
44: MatStructure matstruc;
45: PetscInt n; /* Number of columns in the matrix */
46: PetscInt *colptr; /* Index of first element of each column in row and val */
47: PetscInt *row; /* Row of each element of the matrix */
48: PetscScalar *val; /* Value of each element of the matrix */
49: PetscInt *perm; /* Permutation tabular */
50: PetscInt *invp; /* Reverse permutation tabular */
51: PetscScalar *rhs; /* Rhight-hand-side member */
52: PetscInt rhsnbr; /* Rhight-hand-side number (must be 1) */
53: PetscInt iparm[64]; /* Integer parameters */
54: double dparm[64]; /* Floating point parameters */
55: MPI_Comm pastix_comm; /* PaStiX MPI communicator */
56: PetscMPIInt commRank; /* MPI rank */
57: PetscMPIInt commSize; /* MPI communicator size */
58: PetscBool CleanUpPastix; /* Boolean indicating if we call PaStiX clean step */
59: VecScatter scat_rhs;
60: VecScatter scat_sol;
61: Vec b_seq;
62: PetscBool isAIJ;
63: PetscErrorCode (*Destroy)(Mat);
64: } Mat_Pastix;
66: extern PetscErrorCode MatDuplicate_Pastix(Mat,MatDuplicateOption,Mat*);
70: /*
71: convert Petsc seqaij matrix to CSC: colptr[n], row[nz], val[nz]
73: input:
74: A - matrix in seqaij or mpisbaij (bs=1) format
75: valOnly - FALSE: spaces are allocated and values are set for the CSC
76: TRUE: Only fill values
77: output:
78: n - Size of the matrix
79: colptr - Index of first element of each column in row and val
80: row - Row of each element of the matrix
81: values - Value of each element of the matrix
82: */
83: PetscErrorCode MatConvertToCSC(Mat A,PetscBool valOnly,PetscInt *n,PetscInt **colptr,PetscInt **row,PetscScalar **values)
84: {
85: Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data;
86: PetscInt *rowptr = aa->i;
87: PetscInt *col = aa->j;
88: PetscScalar *rvalues = aa->a;
89: PetscInt m = A->rmap->N;
90: PetscInt nnz;
91: PetscInt i,j, k;
92: PetscInt base = 1;
93: PetscInt idx;
95: PetscInt colidx;
96: PetscInt *colcount;
97: PetscBool isSBAIJ;
98: PetscBool isSeqSBAIJ;
99: PetscBool isMpiSBAIJ;
100: PetscBool isSym;
101: PetscBool flg;
102: PetscInt icntl;
103: PetscInt verb;
104: PetscInt check;
107: MatIsSymmetric(A,0.0,&isSym);
108: PetscObjectTypeCompare((PetscObject)A,MATSBAIJ,&isSBAIJ);
109: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
110: PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMpiSBAIJ);
112: *n = A->cmap->N;
114: /* PaStiX only needs triangular matrix if matrix is symmetric
115: */
116: if (isSym && !(isSBAIJ || isSeqSBAIJ || isMpiSBAIJ)) nnz = (aa->nz - *n)/2 + *n;
117: else nnz = aa->nz;
119: if (!valOnly) {
120: PetscMalloc(((*n)+1) *sizeof(PetscInt) ,colptr);
121: PetscMalloc(nnz *sizeof(PetscInt) ,row);
122: PetscMalloc(nnz *sizeof(PetscScalar),values);
124: if (isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) {
125: PetscMemcpy (*colptr, rowptr, ((*n)+1)*sizeof(PetscInt));
126: for (i = 0; i < *n+1; i++) (*colptr)[i] += base;
127: PetscMemcpy (*row, col, (nnz)*sizeof(PetscInt));
128: for (i = 0; i < nnz; i++) (*row)[i] += base;
129: PetscMemcpy (*values, rvalues, (nnz)*sizeof(PetscScalar));
130: } else {
131: PetscMalloc((*n)*sizeof(PetscInt),&colcount);
133: for (i = 0; i < m; i++) colcount[i] = 0;
134: /* Fill-in colptr */
135: for (i = 0; i < m; i++) {
136: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
137: if (!isSym || col[j] <= i) colcount[col[j]]++;
138: }
139: }
141: (*colptr)[0] = base;
142: for (j = 0; j < *n; j++) {
143: (*colptr)[j+1] = (*colptr)[j] + colcount[j];
144: /* in next loop we fill starting from (*colptr)[colidx] - base */
145: colcount[j] = -base;
146: }
148: /* Fill-in rows and values */
149: for (i = 0; i < m; i++) {
150: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
151: if (!isSym || col[j] <= i) {
152: colidx = col[j];
153: idx = (*colptr)[colidx] + colcount[colidx];
154: (*row)[idx] = i + base;
155: (*values)[idx] = rvalues[j];
156: colcount[colidx]++;
157: }
158: }
159: }
160: PetscFree(colcount);
161: }
162: } else {
163: /* Fill-in only values */
164: for (i = 0; i < m; i++) {
165: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
166: colidx = col[j];
167: if ((isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) ||!isSym || col[j] <= i) {
168: /* look for the value to fill */
169: for (k = (*colptr)[colidx] - base; k < (*colptr)[colidx + 1] - base; k++) {
170: if (((*row)[k]-base) == i) {
171: (*values)[k] = rvalues[j];
172: break;
173: }
174: }
175: /* data structure of sparse matrix has changed */
176: if (k == (*colptr)[colidx + 1] - base) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"overflow on k %D",k);
177: }
178: }
179: }
180: }
182: icntl =-1;
183: check = 0;
184: PetscOptionsInt("-mat_pastix_check","Check the matrix 0 : no, 1 : yes)","None",check,&icntl,&flg);
185: if ((flg && icntl >= 0) || PetscLogPrintInfo) check = icntl;
187: if (check == 1) {
188: PetscScalar *tmpvalues;
189: PetscInt *tmprows,*tmpcolptr;
190: tmpvalues = (PetscScalar*)malloc(nnz*sizeof(PetscScalar)); if (!tmpvalues) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");
191: tmprows = (PetscInt*) malloc(nnz*sizeof(PetscInt)); if (!tmprows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");
192: tmpcolptr = (PetscInt*) malloc((*n+1)*sizeof(PetscInt)); if (!tmpcolptr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");
194: PetscMemcpy(tmpcolptr,*colptr,(*n+1)*sizeof(PetscInt));
195: PetscMemcpy(tmprows,*row,nnz*sizeof(PetscInt));
196: PetscMemcpy(tmpvalues,*values,nnz*sizeof(PetscScalar));
197: PetscFree(*row);
198: PetscFree(*values);
200: icntl=-1;
201: verb = API_VERBOSE_NOT;
202: PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",verb,&icntl,&flg);
203: if ((flg && icntl >= 0) || PetscLogPrintInfo) {
204: verb = icntl;
205: }
206: PASTIX_CHECKMATRIX(MPI_COMM_WORLD,verb,((isSym != 0) ? API_SYM_YES : API_SYM_NO),API_YES,*n,&tmpcolptr,&tmprows,(PastixScalar**)&tmpvalues,NULL,1);
208: PetscMemcpy(*colptr,tmpcolptr,(*n+1)*sizeof(PetscInt));
209: PetscMalloc(((*colptr)[*n]-1)*sizeof(PetscInt),row);
210: PetscMemcpy(*row,tmprows,((*colptr)[*n]-1)*sizeof(PetscInt));
211: PetscMalloc(((*colptr)[*n]-1)*sizeof(PetscScalar),values);
212: PetscMemcpy(*values,tmpvalues,((*colptr)[*n]-1)*sizeof(PetscScalar));
213: free(tmpvalues);
214: free(tmprows);
215: free(tmpcolptr);
217: }
218: return(0);
219: }
225: /*
226: Call clean step of PaStiX if lu->CleanUpPastix == true.
227: Free the CSC matrix.
228: */
229: PetscErrorCode MatDestroy_Pastix(Mat A)
230: {
231: Mat_Pastix *lu=(Mat_Pastix*)A->spptr;
233: PetscMPIInt size=lu->commSize;
236: if (lu && lu->CleanUpPastix) {
237: /* Terminate instance, deallocate memories */
238: if (size > 1) {
239: VecScatterDestroy(&lu->scat_rhs);
240: VecDestroy(&lu->b_seq);
241: VecScatterDestroy(&lu->scat_sol);
242: }
244: lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN;
245: lu->iparm[IPARM_END_TASK] =API_TASK_CLEAN;
247: PASTIX_CALL(&(lu->pastix_data),
248: lu->pastix_comm,
249: lu->n,
250: lu->colptr,
251: lu->row,
252: (PastixScalar*)lu->val,
253: lu->perm,
254: lu->invp,
255: (PastixScalar*)lu->rhs,
256: lu->rhsnbr,
257: lu->iparm,
258: lu->dparm);
260: PetscFree(lu->colptr);
261: PetscFree(lu->row);
262: PetscFree(lu->val);
263: PetscFree(lu->perm);
264: PetscFree(lu->invp);
265: MPI_Comm_free(&(lu->pastix_comm));
266: }
267: if (lu && lu->Destroy) {
268: (lu->Destroy)(A);
269: }
270: PetscFree(A->spptr);
271: return(0);
272: }
276: /*
277: Gather right-hand-side.
278: Call for Solve step.
279: Scatter solution.
280: */
281: PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
282: {
283: Mat_Pastix *lu=(Mat_Pastix*)A->spptr;
284: PetscScalar *array;
285: Vec x_seq;
289: lu->rhsnbr = 1;
290: x_seq = lu->b_seq;
291: if (lu->commSize > 1) {
292: /* PaStiX only supports centralized rhs. Scatter b into a seqential rhs vector */
293: VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
294: VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
295: VecGetArray(x_seq,&array);
296: } else { /* size == 1 */
297: VecCopy(b,x);
298: VecGetArray(x,&array);
299: }
300: lu->rhs = array;
301: if (lu->commSize == 1) {
302: VecRestoreArray(x,&array);
303: } else {
304: VecRestoreArray(x_seq,&array);
305: }
307: /* solve phase */
308: /*-------------*/
309: lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
310: lu->iparm[IPARM_END_TASK] = API_TASK_REFINE;
311: lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;
313: PASTIX_CALL(&(lu->pastix_data),
314: lu->pastix_comm,
315: lu->n,
316: lu->colptr,
317: lu->row,
318: (PastixScalar*)lu->val,
319: lu->perm,
320: lu->invp,
321: (PastixScalar*)lu->rhs,
322: lu->rhsnbr,
323: lu->iparm,
324: lu->dparm);
326: if (lu->iparm[IPARM_ERROR_NUMBER] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in solve phase: lu->iparm[IPARM_ERROR_NUMBER] = %d\n",lu->iparm[IPARM_ERROR_NUMBER]);
328: if (lu->commSize == 1) {
329: VecRestoreArray(x,&(lu->rhs));
330: } else {
331: VecRestoreArray(x_seq,&(lu->rhs));
332: }
334: if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
335: VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
336: VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
337: }
338: return(0);
339: }
341: /*
342: Numeric factorisation using PaStiX solver.
344: */
347: PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
348: {
349: Mat_Pastix *lu =(Mat_Pastix*)(F)->spptr;
350: Mat *tseq;
351: PetscErrorCode 0;
352: PetscInt icntl;
353: PetscInt M=A->rmap->N;
354: PetscBool valOnly,flg, isSym;
355: Mat F_diag;
356: IS is_iden;
357: Vec b;
358: IS isrow;
359: PetscBool isSeqAIJ,isSeqSBAIJ,isMPIAIJ;
362: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
363: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
364: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
365: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
366: (F)->ops->solve = MatSolve_PaStiX;
368: /* Initialize a PASTIX instance */
369: MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm));
370: MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
371: MPI_Comm_size(lu->pastix_comm, &lu->commSize);
373: /* Set pastix options */
374: lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
375: lu->iparm[IPARM_START_TASK] = API_TASK_INIT;
376: lu->iparm[IPARM_END_TASK] = API_TASK_INIT;
378: lu->rhsnbr = 1;
380: /* Call to set default pastix options */
381: PASTIX_CALL(&(lu->pastix_data),
382: lu->pastix_comm,
383: lu->n,
384: lu->colptr,
385: lu->row,
386: (PastixScalar*)lu->val,
387: lu->perm,
388: lu->invp,
389: (PastixScalar*)lu->rhs,
390: lu->rhsnbr,
391: lu->iparm,
392: lu->dparm);
394: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"PaStiX Options","Mat");
396: icntl = -1;
398: lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
400: PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg);
401: if ((flg && icntl >= 0) || PetscLogPrintInfo) {
402: lu->iparm[IPARM_VERBOSE] = icntl;
403: }
404: icntl=-1;
405: PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg);
406: if ((flg && icntl > 0)) {
407: lu->iparm[IPARM_THREAD_NBR] = icntl;
408: }
409: PetscOptionsEnd();
410: valOnly = PETSC_FALSE;
411: } else {
412: if (isSeqAIJ || isMPIAIJ) {
413: PetscFree(lu->colptr);
414: PetscFree(lu->row);
415: PetscFree(lu->val);
416: valOnly = PETSC_FALSE;
417: } else valOnly = PETSC_TRUE;
418: }
420: lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES;
422: /* convert mpi A to seq mat A */
423: ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);
424: MatGetSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);
425: ISDestroy(&isrow);
427: MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
428: MatIsSymmetric(*tseq,0.0,&isSym);
429: MatDestroyMatrices(1,&tseq);
431: if (!lu->perm) {
432: PetscMalloc((lu->n)*sizeof(PetscInt),&(lu->perm));
433: PetscMalloc((lu->n)*sizeof(PetscInt),&(lu->invp));
434: }
436: if (isSym) {
437: /* On symmetric matrix, LLT */
438: lu->iparm[IPARM_SYM] = API_SYM_YES;
439: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
440: } else {
441: /* On unsymmetric matrix, LU */
442: lu->iparm[IPARM_SYM] = API_SYM_NO;
443: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
444: }
446: /*----------------*/
447: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
448: if (!(isSeqAIJ || isSeqSBAIJ)) {
449: /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
450: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);
451: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
452: VecCreate(PetscObjectComm((PetscObject)A),&b);
453: VecSetSizes(b,A->rmap->n,PETSC_DECIDE);
454: VecSetFromOptions(b);
456: VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);
457: VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol);
458: ISDestroy(&is_iden);
459: VecDestroy(&b);
460: }
461: lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
462: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
464: PASTIX_CALL(&(lu->pastix_data),
465: lu->pastix_comm,
466: lu->n,
467: lu->colptr,
468: lu->row,
469: (PastixScalar*)lu->val,
470: lu->perm,
471: lu->invp,
472: (PastixScalar*)lu->rhs,
473: lu->rhsnbr,
474: lu->iparm,
475: lu->dparm);
476: if (lu->iparm[IPARM_ERROR_NUMBER] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in analysis phase: iparm(IPARM_ERROR_NUMBER)=%d\n",lu->iparm[IPARM_ERROR_NUMBER]);
477: } else {
478: lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
479: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
480: PASTIX_CALL(&(lu->pastix_data),
481: lu->pastix_comm,
482: lu->n,
483: lu->colptr,
484: lu->row,
485: (PastixScalar*)lu->val,
486: lu->perm,
487: lu->invp,
488: (PastixScalar*)lu->rhs,
489: lu->rhsnbr,
490: lu->iparm,
491: lu->dparm);
493: if (lu->iparm[IPARM_ERROR_NUMBER] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in analysis phase: iparm(IPARM_ERROR_NUMBER)=%d\n",lu->iparm[IPARM_ERROR_NUMBER]);
494: }
496: if (lu->commSize > 1) {
497: if ((F)->factortype == MAT_FACTOR_LU) {
498: F_diag = ((Mat_MPIAIJ*)(F)->data)->A;
499: } else {
500: F_diag = ((Mat_MPISBAIJ*)(F)->data)->A;
501: }
502: F_diag->assembled = PETSC_TRUE;
503: }
504: (F)->assembled = PETSC_TRUE;
505: lu->matstruc = SAME_NONZERO_PATTERN;
506: lu->CleanUpPastix = PETSC_TRUE;
507: return(0);
508: }
510: /* Note the Petsc r and c permutations are ignored */
513: PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
514: {
515: Mat_Pastix *lu = (Mat_Pastix*)F->spptr;
518: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
519: lu->iparm[IPARM_SYM] = API_SYM_YES;
520: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
521: F->ops->lufactornumeric = MatFactorNumeric_PaStiX;
522: return(0);
523: }
526: /* Note the Petsc r permutation is ignored */
529: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
530: {
531: Mat_Pastix *lu = (Mat_Pastix*)(F)->spptr;
534: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT;
535: lu->iparm[IPARM_SYM] = API_SYM_NO;
536: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
537: (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
538: return(0);
539: }
543: PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
544: {
545: PetscErrorCode ierr;
546: PetscBool iascii;
547: PetscViewerFormat format;
550: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
551: if (iascii) {
552: PetscViewerGetFormat(viewer,&format);
553: if (format == PETSC_VIEWER_ASCII_INFO) {
554: Mat_Pastix *lu=(Mat_Pastix*)A->spptr;
556: PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");
557: PetscViewerASCIIPrintf(viewer," Matrix type : %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
558: PetscViewerASCIIPrintf(viewer," Level of printing (0,1,2): %d \n",lu->iparm[IPARM_VERBOSE]);
559: PetscViewerASCIIPrintf(viewer," Number of refinements iterations : %d \n",lu->iparm[IPARM_NBITER]);
560: PetscPrintf(PETSC_COMM_SELF," Error : %g \n",lu->dparm[DPARM_RELATIVE_ERROR]);
561: }
562: }
563: return(0);
564: }
567: /*MC
568: MATSOLVERPASTIX - A solver package providing direct solvers (LU) for distributed
569: and sequential matrices via the external package PaStiX.
571: Use ./configure --download-pastix to have PETSc installed with PaStiX
573: Options Database Keys:
574: + -mat_pastix_verbose <0,1,2> - print level
575: - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task.
577: Level: beginner
579: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
581: M*/
586: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
587: {
588: Mat_Pastix *lu =(Mat_Pastix*)A->spptr;
591: info->block_size = 1.0;
592: info->nz_allocated = lu->iparm[IPARM_NNZEROS];
593: info->nz_used = lu->iparm[IPARM_NNZEROS];
594: info->nz_unneeded = 0.0;
595: info->assemblies = 0.0;
596: info->mallocs = 0.0;
597: info->memory = 0.0;
598: info->fill_ratio_given = 0;
599: info->fill_ratio_needed = 0;
600: info->factor_mallocs = 0;
601: return(0);
602: }
606: PetscErrorCode MatFactorGetSolverPackage_pastix(Mat A,const MatSolverPackage *type)
607: {
609: *type = MATSOLVERPASTIX;
610: return(0);
611: }
613: /*
614: The seq and mpi versions of this function are the same
615: */
618: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
619: {
620: Mat B;
622: Mat_Pastix *pastix;
625: if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
626: /* Create the factorization matrix */
627: MatCreate(PetscObjectComm((PetscObject)A),&B);
628: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
629: MatSetType(B,((PetscObject)A)->type_name);
630: MatSeqAIJSetPreallocation(B,0,NULL);
632: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
633: B->ops->view = MatView_PaStiX;
634: B->ops->getinfo = MatGetInfo_PaStiX;
636: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
638: B->factortype = MAT_FACTOR_LU;
640: PetscNewLog(B,Mat_Pastix,&pastix);
642: pastix->CleanUpPastix = PETSC_FALSE;
643: pastix->isAIJ = PETSC_TRUE;
644: pastix->scat_rhs = NULL;
645: pastix->scat_sol = NULL;
646: pastix->Destroy = B->ops->destroy;
647: B->ops->destroy = MatDestroy_Pastix;
648: B->spptr = (void*)pastix;
650: *F = B;
651: return(0);
652: }
656: PETSC_EXTERN PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
657: {
658: Mat B;
660: Mat_Pastix *pastix;
663: if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
664: /* Create the factorization matrix */
665: MatCreate(PetscObjectComm((PetscObject)A),&B);
666: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
667: MatSetType(B,((PetscObject)A)->type_name);
668: MatSeqAIJSetPreallocation(B,0,NULL);
669: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
671: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
672: B->ops->view = MatView_PaStiX;
674: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
676: B->factortype = MAT_FACTOR_LU;
678: PetscNewLog(B,Mat_Pastix,&pastix);
680: pastix->CleanUpPastix = PETSC_FALSE;
681: pastix->isAIJ = PETSC_TRUE;
682: pastix->scat_rhs = NULL;
683: pastix->scat_sol = NULL;
684: pastix->Destroy = B->ops->destroy;
685: B->ops->destroy = MatDestroy_Pastix;
686: B->spptr = (void*)pastix;
688: *F = B;
689: return(0);
690: }
694: PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
695: {
696: Mat B;
698: Mat_Pastix *pastix;
701: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
702: /* Create the factorization matrix */
703: MatCreate(PetscObjectComm((PetscObject)A),&B);
704: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
705: MatSetType(B,((PetscObject)A)->type_name);
706: MatSeqSBAIJSetPreallocation(B,1,0,NULL);
707: MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);
709: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
710: B->ops->view = MatView_PaStiX;
712: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
714: B->factortype = MAT_FACTOR_CHOLESKY;
716: PetscNewLog(B,Mat_Pastix,&pastix);
718: pastix->CleanUpPastix = PETSC_FALSE;
719: pastix->isAIJ = PETSC_TRUE;
720: pastix->scat_rhs = NULL;
721: pastix->scat_sol = NULL;
722: pastix->Destroy = B->ops->destroy;
723: B->ops->destroy = MatDestroy_Pastix;
724: B->spptr = (void*)pastix;
726: *F = B;
727: return(0);
728: }
732: PETSC_EXTERN PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
733: {
734: Mat B;
736: Mat_Pastix *pastix;
739: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
741: /* Create the factorization matrix */
742: MatCreate(PetscObjectComm((PetscObject)A),&B);
743: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
744: MatSetType(B,((PetscObject)A)->type_name);
745: MatSeqSBAIJSetPreallocation(B,1,0,NULL);
746: MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);
748: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
749: B->ops->view = MatView_PaStiX;
751: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
753: B->factortype = MAT_FACTOR_CHOLESKY;
755: PetscNewLog(B,Mat_Pastix,&pastix);
757: pastix->CleanUpPastix = PETSC_FALSE;
758: pastix->isAIJ = PETSC_TRUE;
759: pastix->scat_rhs = NULL;
760: pastix->scat_sol = NULL;
761: pastix->Destroy = B->ops->destroy;
762: B->ops->destroy = MatDestroy_Pastix;
763: B->spptr = (void*)pastix;
765: *F = B;
766: return(0);
767: }