Actual source code: pastix.c
petsc-3.5.4 2015-05-23
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)
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: #else
22: #define PASTIX_CALL z_pastix
23: #define PASTIX_CHECKMATRIX z_pastix_checkMatrix
24: #endif
26: #else /* PETSC_USE_COMPLEX */
28: #if defined(PETSC_USE_REAL_SINGLE)
29: #define PASTIX_CALL s_pastix
30: #define PASTIX_CHECKMATRIX s_pastix_checkMatrix
31: #else
32: #define PASTIX_CALL d_pastix
33: #define PASTIX_CHECKMATRIX d_pastix_checkMatrix
34: #endif
36: #endif /* PETSC_USE_COMPLEX */
38: typedef PetscScalar PastixScalar;
40: typedef struct Mat_Pastix_ {
41: pastix_data_t *pastix_data; /* Pastix data storage structure */
42: MatStructure matstruc;
43: PetscInt n; /* Number of columns in the matrix */
44: PetscInt *colptr; /* Index of first element of each column in row and val */
45: PetscInt *row; /* Row of each element of the matrix */
46: PetscScalar *val; /* Value of each element of the matrix */
47: PetscInt *perm; /* Permutation tabular */
48: PetscInt *invp; /* Reverse permutation tabular */
49: PetscScalar *rhs; /* Rhight-hand-side member */
50: PetscInt rhsnbr; /* Rhight-hand-side number (must be 1) */
51: PetscInt iparm[64]; /* Integer parameters */
52: double dparm[64]; /* Floating point parameters */
53: MPI_Comm pastix_comm; /* PaStiX MPI communicator */
54: PetscMPIInt commRank; /* MPI rank */
55: PetscMPIInt commSize; /* MPI communicator size */
56: PetscBool CleanUpPastix; /* Boolean indicating if we call PaStiX clean step */
57: VecScatter scat_rhs;
58: VecScatter scat_sol;
59: Vec b_seq;
60: PetscBool isAIJ;
61: PetscErrorCode (*Destroy)(Mat);
62: } Mat_Pastix;
64: extern PetscErrorCode MatDuplicate_Pastix(Mat,MatDuplicateOption,Mat*);
68: /*
69: convert Petsc seqaij matrix to CSC: colptr[n], row[nz], val[nz]
71: input:
72: A - matrix in seqaij or mpisbaij (bs=1) format
73: valOnly - FALSE: spaces are allocated and values are set for the CSC
74: TRUE: Only fill values
75: output:
76: n - Size of the matrix
77: colptr - Index of first element of each column in row and val
78: row - Row of each element of the matrix
79: values - Value of each element of the matrix
80: */
81: PetscErrorCode MatConvertToCSC(Mat A,PetscBool valOnly,PetscInt *n,PetscInt **colptr,PetscInt **row,PetscScalar **values)
82: {
83: Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data;
84: PetscInt *rowptr = aa->i;
85: PetscInt *col = aa->j;
86: PetscScalar *rvalues = aa->a;
87: PetscInt m = A->rmap->N;
88: PetscInt nnz;
89: PetscInt i,j, k;
90: PetscInt base = 1;
91: PetscInt idx;
93: PetscInt colidx;
94: PetscInt *colcount;
95: PetscBool isSBAIJ;
96: PetscBool isSeqSBAIJ;
97: PetscBool isMpiSBAIJ;
98: PetscBool isSym;
99: PetscBool flg;
100: PetscInt icntl;
101: PetscInt verb;
102: PetscInt check;
105: MatIsSymmetric(A,0.0,&isSym);
106: PetscObjectTypeCompare((PetscObject)A,MATSBAIJ,&isSBAIJ);
107: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
108: PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMpiSBAIJ);
110: *n = A->cmap->N;
112: /* PaStiX only needs triangular matrix if matrix is symmetric
113: */
114: if (isSym && !(isSBAIJ || isSeqSBAIJ || isMpiSBAIJ)) nnz = (aa->nz - *n)/2 + *n;
115: else nnz = aa->nz;
117: if (!valOnly) {
118: PetscMalloc(((*n)+1) *sizeof(PetscInt) ,colptr);
119: PetscMalloc(nnz *sizeof(PetscInt) ,row);
120: PetscMalloc1(nnz ,values);
122: if (isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) {
123: PetscMemcpy (*colptr, rowptr, ((*n)+1)*sizeof(PetscInt));
124: for (i = 0; i < *n+1; i++) (*colptr)[i] += base;
125: PetscMemcpy (*row, col, (nnz)*sizeof(PetscInt));
126: for (i = 0; i < nnz; i++) (*row)[i] += base;
127: PetscMemcpy (*values, rvalues, (nnz)*sizeof(PetscScalar));
128: } else {
129: PetscMalloc1((*n),&colcount);
131: for (i = 0; i < m; i++) colcount[i] = 0;
132: /* Fill-in colptr */
133: for (i = 0; i < m; i++) {
134: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
135: if (!isSym || col[j] <= i) colcount[col[j]]++;
136: }
137: }
139: (*colptr)[0] = base;
140: for (j = 0; j < *n; j++) {
141: (*colptr)[j+1] = (*colptr)[j] + colcount[j];
142: /* in next loop we fill starting from (*colptr)[colidx] - base */
143: colcount[j] = -base;
144: }
146: /* Fill-in rows and values */
147: for (i = 0; i < m; i++) {
148: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
149: if (!isSym || col[j] <= i) {
150: colidx = col[j];
151: idx = (*colptr)[colidx] + colcount[colidx];
152: (*row)[idx] = i + base;
153: (*values)[idx] = rvalues[j];
154: colcount[colidx]++;
155: }
156: }
157: }
158: PetscFree(colcount);
159: }
160: } else {
161: /* Fill-in only values */
162: for (i = 0; i < m; i++) {
163: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
164: colidx = col[j];
165: if ((isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) ||!isSym || col[j] <= i) {
166: /* look for the value to fill */
167: for (k = (*colptr)[colidx] - base; k < (*colptr)[colidx + 1] - base; k++) {
168: if (((*row)[k]-base) == i) {
169: (*values)[k] = rvalues[j];
170: break;
171: }
172: }
173: /* data structure of sparse matrix has changed */
174: if (k == (*colptr)[colidx + 1] - base) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"overflow on k %D",k);
175: }
176: }
177: }
178: }
180: icntl =-1;
181: check = 0;
182: PetscOptionsGetInt(((PetscObject) A)->prefix, "-mat_pastix_check", &icntl, &flg);
183: if ((flg && icntl >= 0) || PetscLogPrintInfo) check = icntl;
185: if (check == 1) {
186: PetscScalar *tmpvalues;
187: PetscInt *tmprows,*tmpcolptr;
188: tmpvalues = (PetscScalar*)malloc(nnz*sizeof(PetscScalar)); if (!tmpvalues) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");
189: tmprows = (PetscInt*) malloc(nnz*sizeof(PetscInt)); if (!tmprows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");
190: tmpcolptr = (PetscInt*) malloc((*n+1)*sizeof(PetscInt)); if (!tmpcolptr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");
192: PetscMemcpy(tmpcolptr,*colptr,(*n+1)*sizeof(PetscInt));
193: PetscMemcpy(tmprows,*row,nnz*sizeof(PetscInt));
194: PetscMemcpy(tmpvalues,*values,nnz*sizeof(PetscScalar));
195: PetscFree(*row);
196: PetscFree(*values);
198: icntl=-1;
199: verb = API_VERBOSE_NOT;
200: /* "iparm[IPARM_VERBOSE] : level of printing (0 to 2)" */
201: PetscOptionsGetInt(((PetscObject) A)->prefix, "-mat_pastix_verbose", &icntl, &flg);
202: if ((flg && icntl >= 0) || PetscLogPrintInfo) verb = icntl;
203: PASTIX_CHECKMATRIX(MPI_COMM_WORLD,verb,((isSym != 0) ? API_SYM_YES : API_SYM_NO),API_YES,*n,&tmpcolptr,&tmprows,(PastixScalar**)&tmpvalues,NULL,1);
205: PetscMemcpy(*colptr,tmpcolptr,(*n+1)*sizeof(PetscInt));
206: PetscMalloc1(((*colptr)[*n]-1),row);
207: PetscMemcpy(*row,tmprows,((*colptr)[*n]-1)*sizeof(PetscInt));
208: PetscMalloc1(((*colptr)[*n]-1),values);
209: PetscMemcpy(*values,tmpvalues,((*colptr)[*n]-1)*sizeof(PetscScalar));
210: free(tmpvalues);
211: free(tmprows);
212: free(tmpcolptr);
214: }
215: return(0);
216: }
222: /*
223: Call clean step of PaStiX if lu->CleanUpPastix == true.
224: Free the CSC matrix.
225: */
226: PetscErrorCode MatDestroy_Pastix(Mat A)
227: {
228: Mat_Pastix *lu=(Mat_Pastix*)A->spptr;
230: PetscMPIInt size=lu->commSize;
233: if (lu && lu->CleanUpPastix) {
234: /* Terminate instance, deallocate memories */
235: if (size > 1) {
236: VecScatterDestroy(&lu->scat_rhs);
237: VecDestroy(&lu->b_seq);
238: VecScatterDestroy(&lu->scat_sol);
239: }
241: lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN;
242: lu->iparm[IPARM_END_TASK] =API_TASK_CLEAN;
244: PASTIX_CALL(&(lu->pastix_data),
245: lu->pastix_comm,
246: lu->n,
247: lu->colptr,
248: lu->row,
249: (PastixScalar*)lu->val,
250: lu->perm,
251: lu->invp,
252: (PastixScalar*)lu->rhs,
253: lu->rhsnbr,
254: lu->iparm,
255: lu->dparm);
257: PetscFree(lu->colptr);
258: PetscFree(lu->row);
259: PetscFree(lu->val);
260: PetscFree(lu->perm);
261: PetscFree(lu->invp);
262: MPI_Comm_free(&(lu->pastix_comm));
263: }
264: if (lu && lu->Destroy) {
265: (lu->Destroy)(A);
266: }
267: PetscFree(A->spptr);
268: return(0);
269: }
273: /*
274: Gather right-hand-side.
275: Call for Solve step.
276: Scatter solution.
277: */
278: PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
279: {
280: Mat_Pastix *lu=(Mat_Pastix*)A->spptr;
281: PetscScalar *array;
282: Vec x_seq;
286: lu->rhsnbr = 1;
287: x_seq = lu->b_seq;
288: if (lu->commSize > 1) {
289: /* PaStiX only supports centralized rhs. Scatter b into a seqential rhs vector */
290: VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
291: VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
292: VecGetArray(x_seq,&array);
293: } else { /* size == 1 */
294: VecCopy(b,x);
295: VecGetArray(x,&array);
296: }
297: lu->rhs = array;
298: if (lu->commSize == 1) {
299: VecRestoreArray(x,&array);
300: } else {
301: VecRestoreArray(x_seq,&array);
302: }
304: /* solve phase */
305: /*-------------*/
306: lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
307: lu->iparm[IPARM_END_TASK] = API_TASK_REFINE;
308: lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;
310: PASTIX_CALL(&(lu->pastix_data),
311: lu->pastix_comm,
312: lu->n,
313: lu->colptr,
314: lu->row,
315: (PastixScalar*)lu->val,
316: lu->perm,
317: lu->invp,
318: (PastixScalar*)lu->rhs,
319: lu->rhsnbr,
320: lu->iparm,
321: lu->dparm);
323: 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]);
325: if (lu->commSize == 1) {
326: VecRestoreArray(x,&(lu->rhs));
327: } else {
328: VecRestoreArray(x_seq,&(lu->rhs));
329: }
331: if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
332: VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
333: VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
334: }
335: return(0);
336: }
338: /*
339: Numeric factorisation using PaStiX solver.
341: */
344: PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
345: {
346: Mat_Pastix *lu =(Mat_Pastix*)(F)->spptr;
347: Mat *tseq;
348: PetscErrorCode 0;
349: PetscInt icntl;
350: PetscInt M=A->rmap->N;
351: PetscBool valOnly,flg, isSym;
352: Mat F_diag;
353: IS is_iden;
354: Vec b;
355: IS isrow;
356: PetscBool isSeqAIJ,isSeqSBAIJ,isMPIAIJ;
359: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
360: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
361: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
362: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
363: (F)->ops->solve = MatSolve_PaStiX;
365: /* Initialize a PASTIX instance */
366: MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm));
367: MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
368: MPI_Comm_size(lu->pastix_comm, &lu->commSize);
370: /* Set pastix options */
371: lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
372: lu->iparm[IPARM_START_TASK] = API_TASK_INIT;
373: lu->iparm[IPARM_END_TASK] = API_TASK_INIT;
375: lu->rhsnbr = 1;
377: /* Call to set default pastix options */
378: PASTIX_CALL(&(lu->pastix_data),
379: lu->pastix_comm,
380: lu->n,
381: lu->colptr,
382: lu->row,
383: (PastixScalar*)lu->val,
384: lu->perm,
385: lu->invp,
386: (PastixScalar*)lu->rhs,
387: lu->rhsnbr,
388: lu->iparm,
389: lu->dparm);
391: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"PaStiX Options","Mat");
393: icntl = -1;
395: lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
397: PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg);
398: if ((flg && icntl >= 0) || PetscLogPrintInfo) {
399: lu->iparm[IPARM_VERBOSE] = icntl;
400: }
401: icntl=-1;
402: PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg);
403: if ((flg && icntl > 0)) {
404: lu->iparm[IPARM_THREAD_NBR] = icntl;
405: }
406: PetscOptionsEnd();
407: valOnly = PETSC_FALSE;
408: } else {
409: if (isSeqAIJ || isMPIAIJ) {
410: PetscFree(lu->colptr);
411: PetscFree(lu->row);
412: PetscFree(lu->val);
413: valOnly = PETSC_FALSE;
414: } else valOnly = PETSC_TRUE;
415: }
417: lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES;
419: /* convert mpi A to seq mat A */
420: ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);
421: MatGetSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);
422: ISDestroy(&isrow);
424: MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
425: MatIsSymmetric(*tseq,0.0,&isSym);
426: MatDestroyMatrices(1,&tseq);
428: if (!lu->perm) {
429: PetscMalloc1((lu->n),&(lu->perm));
430: PetscMalloc1((lu->n),&(lu->invp));
431: }
433: if (isSym) {
434: /* On symmetric matrix, LLT */
435: lu->iparm[IPARM_SYM] = API_SYM_YES;
436: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
437: } else {
438: /* On unsymmetric matrix, LU */
439: lu->iparm[IPARM_SYM] = API_SYM_NO;
440: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
441: }
443: /*----------------*/
444: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
445: if (!(isSeqAIJ || isSeqSBAIJ)) {
446: /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
447: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);
448: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
449: MatGetVecs(A,NULL,&b);
450: VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);
451: VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol);
452: ISDestroy(&is_iden);
453: VecDestroy(&b);
454: }
455: lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
456: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
458: PASTIX_CALL(&(lu->pastix_data),
459: lu->pastix_comm,
460: lu->n,
461: lu->colptr,
462: lu->row,
463: (PastixScalar*)lu->val,
464: lu->perm,
465: lu->invp,
466: (PastixScalar*)lu->rhs,
467: lu->rhsnbr,
468: lu->iparm,
469: lu->dparm);
470: 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]);
471: } else {
472: lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
473: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
474: PASTIX_CALL(&(lu->pastix_data),
475: lu->pastix_comm,
476: lu->n,
477: lu->colptr,
478: lu->row,
479: (PastixScalar*)lu->val,
480: lu->perm,
481: lu->invp,
482: (PastixScalar*)lu->rhs,
483: lu->rhsnbr,
484: lu->iparm,
485: lu->dparm);
487: 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]);
488: }
490: if (lu->commSize > 1) {
491: if ((F)->factortype == MAT_FACTOR_LU) {
492: F_diag = ((Mat_MPIAIJ*)(F)->data)->A;
493: } else {
494: F_diag = ((Mat_MPISBAIJ*)(F)->data)->A;
495: }
496: F_diag->assembled = PETSC_TRUE;
497: }
498: (F)->assembled = PETSC_TRUE;
499: lu->matstruc = SAME_NONZERO_PATTERN;
500: lu->CleanUpPastix = PETSC_TRUE;
501: return(0);
502: }
504: /* Note the Petsc r and c permutations are ignored */
507: PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
508: {
509: Mat_Pastix *lu = (Mat_Pastix*)F->spptr;
512: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
513: lu->iparm[IPARM_SYM] = API_SYM_YES;
514: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
515: F->ops->lufactornumeric = MatFactorNumeric_PaStiX;
516: return(0);
517: }
520: /* Note the Petsc r permutation is ignored */
523: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
524: {
525: Mat_Pastix *lu = (Mat_Pastix*)(F)->spptr;
528: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT;
529: lu->iparm[IPARM_SYM] = API_SYM_NO;
530: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
531: (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
532: return(0);
533: }
537: PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
538: {
539: PetscErrorCode ierr;
540: PetscBool iascii;
541: PetscViewerFormat format;
544: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
545: if (iascii) {
546: PetscViewerGetFormat(viewer,&format);
547: if (format == PETSC_VIEWER_ASCII_INFO) {
548: Mat_Pastix *lu=(Mat_Pastix*)A->spptr;
550: PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");
551: PetscViewerASCIIPrintf(viewer," Matrix type : %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
552: PetscViewerASCIIPrintf(viewer," Level of printing (0,1,2): %d \n",lu->iparm[IPARM_VERBOSE]);
553: PetscViewerASCIIPrintf(viewer," Number of refinements iterations : %d \n",lu->iparm[IPARM_NBITER]);
554: PetscPrintf(PETSC_COMM_SELF," Error : %g \n",lu->dparm[DPARM_RELATIVE_ERROR]);
555: }
556: }
557: return(0);
558: }
561: /*MC
562: MATSOLVERPASTIX - A solver package providing direct solvers (LU) for distributed
563: and sequential matrices via the external package PaStiX.
565: Use ./configure --download-pastix to have PETSc installed with PaStiX
567: Options Database Keys:
568: + -mat_pastix_verbose <0,1,2> - print level
569: - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task.
571: Level: beginner
573: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
575: M*/
580: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
581: {
582: Mat_Pastix *lu =(Mat_Pastix*)A->spptr;
585: info->block_size = 1.0;
586: info->nz_allocated = lu->iparm[IPARM_NNZEROS];
587: info->nz_used = lu->iparm[IPARM_NNZEROS];
588: info->nz_unneeded = 0.0;
589: info->assemblies = 0.0;
590: info->mallocs = 0.0;
591: info->memory = 0.0;
592: info->fill_ratio_given = 0;
593: info->fill_ratio_needed = 0;
594: info->factor_mallocs = 0;
595: return(0);
596: }
600: PetscErrorCode MatFactorGetSolverPackage_pastix(Mat A,const MatSolverPackage *type)
601: {
603: *type = MATSOLVERPASTIX;
604: return(0);
605: }
607: /*
608: The seq and mpi versions of this function are the same
609: */
612: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
613: {
614: Mat B;
616: Mat_Pastix *pastix;
619: if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
620: /* Create the factorization matrix */
621: MatCreate(PetscObjectComm((PetscObject)A),&B);
622: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
623: MatSetType(B,((PetscObject)A)->type_name);
624: MatSeqAIJSetPreallocation(B,0,NULL);
626: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
627: B->ops->view = MatView_PaStiX;
628: B->ops->getinfo = MatGetInfo_PaStiX;
630: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
632: B->factortype = MAT_FACTOR_LU;
634: PetscNewLog(B,&pastix);
636: pastix->CleanUpPastix = PETSC_FALSE;
637: pastix->isAIJ = PETSC_TRUE;
638: pastix->scat_rhs = NULL;
639: pastix->scat_sol = NULL;
640: pastix->Destroy = B->ops->destroy;
641: B->ops->destroy = MatDestroy_Pastix;
642: B->spptr = (void*)pastix;
644: *F = B;
645: return(0);
646: }
650: PETSC_EXTERN PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
651: {
652: Mat B;
654: Mat_Pastix *pastix;
657: if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
658: /* Create the factorization matrix */
659: MatCreate(PetscObjectComm((PetscObject)A),&B);
660: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
661: MatSetType(B,((PetscObject)A)->type_name);
662: MatSeqAIJSetPreallocation(B,0,NULL);
663: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
665: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
666: B->ops->view = MatView_PaStiX;
668: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
670: B->factortype = MAT_FACTOR_LU;
672: PetscNewLog(B,&pastix);
674: pastix->CleanUpPastix = PETSC_FALSE;
675: pastix->isAIJ = PETSC_TRUE;
676: pastix->scat_rhs = NULL;
677: pastix->scat_sol = NULL;
678: pastix->Destroy = B->ops->destroy;
679: B->ops->destroy = MatDestroy_Pastix;
680: B->spptr = (void*)pastix;
682: *F = B;
683: return(0);
684: }
688: PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
689: {
690: Mat B;
692: Mat_Pastix *pastix;
695: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
696: /* Create the factorization matrix */
697: MatCreate(PetscObjectComm((PetscObject)A),&B);
698: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
699: MatSetType(B,((PetscObject)A)->type_name);
700: MatSeqSBAIJSetPreallocation(B,1,0,NULL);
701: MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);
703: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
704: B->ops->view = MatView_PaStiX;
706: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
708: B->factortype = MAT_FACTOR_CHOLESKY;
710: PetscNewLog(B,&pastix);
712: pastix->CleanUpPastix = PETSC_FALSE;
713: pastix->isAIJ = PETSC_TRUE;
714: pastix->scat_rhs = NULL;
715: pastix->scat_sol = NULL;
716: pastix->Destroy = B->ops->destroy;
717: B->ops->destroy = MatDestroy_Pastix;
718: B->spptr = (void*)pastix;
720: *F = B;
721: return(0);
722: }
726: PETSC_EXTERN PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
727: {
728: Mat B;
730: Mat_Pastix *pastix;
733: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
735: /* Create the factorization matrix */
736: MatCreate(PetscObjectComm((PetscObject)A),&B);
737: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
738: MatSetType(B,((PetscObject)A)->type_name);
739: MatSeqSBAIJSetPreallocation(B,1,0,NULL);
740: MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);
742: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
743: B->ops->view = MatView_PaStiX;
745: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
747: B->factortype = MAT_FACTOR_CHOLESKY;
749: PetscNewLog(B,&pastix);
751: pastix->CleanUpPastix = PETSC_FALSE;
752: pastix->isAIJ = PETSC_TRUE;
753: pastix->scat_rhs = NULL;
754: pastix->scat_sol = NULL;
755: pastix->Destroy = B->ops->destroy;
756: B->ops->destroy = MatDestroy_Pastix;
757: B->spptr = (void*)pastix;
759: *F = B;
760: return(0);
761: }