Actual source code: pastix.c
petsc-3.8.4 2018-03-24
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: } Mat_Pastix;
62: extern PetscErrorCode MatDuplicate_Pastix(Mat,MatDuplicateOption,Mat*);
64: /*
65: convert Petsc seqaij matrix to CSC: colptr[n], row[nz], val[nz]
67: input:
68: A - matrix in seqaij or mpisbaij (bs=1) format
69: valOnly - FALSE: spaces are allocated and values are set for the CSC
70: TRUE: Only fill values
71: output:
72: n - Size of the matrix
73: colptr - Index of first element of each column in row and val
74: row - Row of each element of the matrix
75: values - Value of each element of the matrix
76: */
77: PetscErrorCode MatConvertToCSC(Mat A,PetscBool valOnly,PetscInt *n,PetscInt **colptr,PetscInt **row,PetscScalar **values)
78: {
79: Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data;
80: PetscInt *rowptr = aa->i;
81: PetscInt *col = aa->j;
82: PetscScalar *rvalues = aa->a;
83: PetscInt m = A->rmap->N;
84: PetscInt nnz;
85: PetscInt i,j, k;
86: PetscInt base = 1;
87: PetscInt idx;
89: PetscInt colidx;
90: PetscInt *colcount;
91: PetscBool isSBAIJ;
92: PetscBool isSeqSBAIJ;
93: PetscBool isMpiSBAIJ;
94: PetscBool isSym;
95: PetscBool flg;
96: PetscInt icntl;
97: PetscInt verb;
98: PetscInt check;
101: MatIsSymmetric(A,0.0,&isSym);
102: PetscObjectTypeCompare((PetscObject)A,MATSBAIJ,&isSBAIJ);
103: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
104: PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMpiSBAIJ);
106: *n = A->cmap->N;
108: /* PaStiX only needs triangular matrix if matrix is symmetric
109: */
110: if (isSym && !(isSBAIJ || isSeqSBAIJ || isMpiSBAIJ)) nnz = (aa->nz - *n)/2 + *n;
111: else nnz = aa->nz;
113: if (!valOnly) {
114: PetscMalloc1((*n)+1,colptr);
115: PetscMalloc1(nnz,row);
116: PetscMalloc1(nnz,values);
118: if (isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) {
119: PetscMemcpy (*colptr, rowptr, ((*n)+1)*sizeof(PetscInt));
120: for (i = 0; i < *n+1; i++) (*colptr)[i] += base;
121: PetscMemcpy (*row, col, (nnz)*sizeof(PetscInt));
122: for (i = 0; i < nnz; i++) (*row)[i] += base;
123: PetscMemcpy (*values, rvalues, (nnz)*sizeof(PetscScalar));
124: } else {
125: PetscMalloc1(*n,&colcount);
127: for (i = 0; i < m; i++) colcount[i] = 0;
128: /* Fill-in colptr */
129: for (i = 0; i < m; i++) {
130: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
131: if (!isSym || col[j] <= i) colcount[col[j]]++;
132: }
133: }
135: (*colptr)[0] = base;
136: for (j = 0; j < *n; j++) {
137: (*colptr)[j+1] = (*colptr)[j] + colcount[j];
138: /* in next loop we fill starting from (*colptr)[colidx] - base */
139: colcount[j] = -base;
140: }
142: /* Fill-in rows and values */
143: for (i = 0; i < m; i++) {
144: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
145: if (!isSym || col[j] <= i) {
146: colidx = col[j];
147: idx = (*colptr)[colidx] + colcount[colidx];
148: (*row)[idx] = i + base;
149: (*values)[idx] = rvalues[j];
150: colcount[colidx]++;
151: }
152: }
153: }
154: PetscFree(colcount);
155: }
156: } else {
157: /* Fill-in only values */
158: for (i = 0; i < m; i++) {
159: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
160: colidx = col[j];
161: if ((isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) ||!isSym || col[j] <= i) {
162: /* look for the value to fill */
163: for (k = (*colptr)[colidx] - base; k < (*colptr)[colidx + 1] - base; k++) {
164: if (((*row)[k]-base) == i) {
165: (*values)[k] = rvalues[j];
166: break;
167: }
168: }
169: /* data structure of sparse matrix has changed */
170: if (k == (*colptr)[colidx + 1] - base) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"overflow on k %D",k);
171: }
172: }
173: }
174: }
176: icntl =-1;
177: check = 0;
178: PetscOptionsGetInt(NULL,((PetscObject) A)->prefix, "-mat_pastix_check", &icntl, &flg);
179: if ((flg && icntl >= 0) || PetscLogPrintInfo) check = icntl;
181: if (check == 1) {
182: PetscScalar *tmpvalues;
183: PetscInt *tmprows,*tmpcolptr;
185: PetscMalloc3(nnz,&tmpvalues,nnz,&tmprows,*n+1,&tmpcolptr);
187: PetscMemcpy(tmpcolptr,*colptr,(*n+1)*sizeof(PetscInt));
188: PetscMemcpy(tmprows,*row,nnz*sizeof(PetscInt));
189: PetscMemcpy(tmpvalues,*values,nnz*sizeof(PetscScalar));
190: PetscFree(*row);
191: PetscFree(*values);
193: icntl=-1;
194: verb = API_VERBOSE_NOT;
195: /* "iparm[IPARM_VERBOSE] : level of printing (0 to 2)" */
196: PetscOptionsGetInt(NULL,((PetscObject) A)->prefix, "-mat_pastix_verbose", &icntl, &flg);
197: if ((flg && icntl >= 0) || PetscLogPrintInfo) verb = icntl;
198: PASTIX_CHECKMATRIX(PetscObjectComm((PetscObject)A),verb,((isSym != 0) ? API_SYM_YES : API_SYM_NO),API_YES,*n,&tmpcolptr,&tmprows,(PastixScalar**)&tmpvalues,NULL,1);
200: PetscMemcpy(*colptr,tmpcolptr,(*n+1)*sizeof(PetscInt));
201: PetscMalloc1(((*colptr)[*n]-1),row);
202: PetscMemcpy(*row,tmprows,((*colptr)[*n]-1)*sizeof(PetscInt));
203: PetscMalloc1(((*colptr)[*n]-1),values);
204: PetscMemcpy(*values,tmpvalues,((*colptr)[*n]-1)*sizeof(PetscScalar));
205: PetscFree3(tmpvalues,tmprows,tmpcolptr);
207: }
208: return(0);
209: }
211: /*
212: Call clean step of PaStiX if lu->CleanUpPastix == true.
213: Free the CSC matrix.
214: */
215: PetscErrorCode MatDestroy_Pastix(Mat A)
216: {
217: Mat_Pastix *lu=(Mat_Pastix*)A->data;
221: if (lu->CleanUpPastix) {
222: /* Terminate instance, deallocate memories */
223: VecScatterDestroy(&lu->scat_rhs);
224: VecDestroy(&lu->b_seq);
225: VecScatterDestroy(&lu->scat_sol);
227: lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN;
228: lu->iparm[IPARM_END_TASK] =API_TASK_CLEAN;
230: PASTIX_CALL(&(lu->pastix_data),
231: lu->pastix_comm,
232: lu->n,
233: lu->colptr,
234: lu->row,
235: (PastixScalar*)lu->val,
236: lu->perm,
237: lu->invp,
238: (PastixScalar*)lu->rhs,
239: lu->rhsnbr,
240: lu->iparm,
241: lu->dparm);
243: PetscFree(lu->colptr);
244: PetscFree(lu->row);
245: PetscFree(lu->val);
246: PetscFree(lu->perm);
247: PetscFree(lu->invp);
248: MPI_Comm_free(&(lu->pastix_comm));
249: }
250: PetscFree(A->data);
251: return(0);
252: }
254: /*
255: Gather right-hand-side.
256: Call for Solve step.
257: Scatter solution.
258: */
259: PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
260: {
261: Mat_Pastix *lu=(Mat_Pastix*)A->data;
262: PetscScalar *array;
263: Vec x_seq;
267: lu->rhsnbr = 1;
268: x_seq = lu->b_seq;
269: if (lu->commSize > 1) {
270: /* PaStiX only supports centralized rhs. Scatter b into a seqential rhs vector */
271: VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
272: VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
273: VecGetArray(x_seq,&array);
274: } else { /* size == 1 */
275: VecCopy(b,x);
276: VecGetArray(x,&array);
277: }
278: lu->rhs = array;
279: if (lu->commSize == 1) {
280: VecRestoreArray(x,&array);
281: } else {
282: VecRestoreArray(x_seq,&array);
283: }
285: /* solve phase */
286: /*-------------*/
287: lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
288: lu->iparm[IPARM_END_TASK] = API_TASK_REFINE;
289: lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;
291: PASTIX_CALL(&(lu->pastix_data),
292: lu->pastix_comm,
293: lu->n,
294: lu->colptr,
295: lu->row,
296: (PastixScalar*)lu->val,
297: lu->perm,
298: lu->invp,
299: (PastixScalar*)lu->rhs,
300: lu->rhsnbr,
301: lu->iparm,
302: lu->dparm);
304: 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]);
306: if (lu->commSize == 1) {
307: VecRestoreArray(x,&(lu->rhs));
308: } else {
309: VecRestoreArray(x_seq,&(lu->rhs));
310: }
312: if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
313: VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
314: VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
315: }
316: return(0);
317: }
319: /*
320: Numeric factorisation using PaStiX solver.
322: */
323: PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
324: {
325: Mat_Pastix *lu =(Mat_Pastix*)(F)->data;
326: Mat *tseq;
327: PetscErrorCode 0;
328: PetscInt icntl;
329: PetscInt M=A->rmap->N;
330: PetscBool valOnly,flg, isSym;
331: IS is_iden;
332: Vec b;
333: IS isrow;
334: PetscBool isSeqAIJ,isSeqSBAIJ,isMPIAIJ;
337: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
338: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
339: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
340: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
341: (F)->ops->solve = MatSolve_PaStiX;
343: /* Initialize a PASTIX instance */
344: MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm));
345: MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
346: MPI_Comm_size(lu->pastix_comm, &lu->commSize);
348: /* Set pastix options */
349: lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
350: lu->iparm[IPARM_START_TASK] = API_TASK_INIT;
351: lu->iparm[IPARM_END_TASK] = API_TASK_INIT;
353: lu->rhsnbr = 1;
355: /* Call to set default pastix options */
356: PASTIX_CALL(&(lu->pastix_data),
357: lu->pastix_comm,
358: lu->n,
359: lu->colptr,
360: lu->row,
361: (PastixScalar*)lu->val,
362: lu->perm,
363: lu->invp,
364: (PastixScalar*)lu->rhs,
365: lu->rhsnbr,
366: lu->iparm,
367: lu->dparm);
369: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"PaStiX Options","Mat");
371: icntl = -1;
373: lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
375: PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg);
376: if ((flg && icntl >= 0) || PetscLogPrintInfo) {
377: lu->iparm[IPARM_VERBOSE] = icntl;
378: }
379: icntl=-1;
380: PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg);
381: if ((flg && icntl > 0)) {
382: lu->iparm[IPARM_THREAD_NBR] = icntl;
383: }
384: PetscOptionsEnd();
385: valOnly = PETSC_FALSE;
386: } else {
387: if (isSeqAIJ || isMPIAIJ) {
388: PetscFree(lu->colptr);
389: PetscFree(lu->row);
390: PetscFree(lu->val);
391: valOnly = PETSC_FALSE;
392: } else valOnly = PETSC_TRUE;
393: }
395: lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES;
397: /* convert mpi A to seq mat A */
398: ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);
399: MatCreateSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);
400: ISDestroy(&isrow);
402: MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
403: MatIsSymmetric(*tseq,0.0,&isSym);
404: MatDestroyMatrices(1,&tseq);
406: if (!lu->perm) {
407: PetscMalloc1(lu->n,&(lu->perm));
408: PetscMalloc1(lu->n,&(lu->invp));
409: }
411: if (isSym) {
412: /* On symmetric matrix, LLT */
413: lu->iparm[IPARM_SYM] = API_SYM_YES;
414: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
415: } else {
416: /* On unsymmetric matrix, LU */
417: lu->iparm[IPARM_SYM] = API_SYM_NO;
418: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
419: }
421: /*----------------*/
422: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
423: if (!(isSeqAIJ || isSeqSBAIJ) && !lu->b_seq) {
424: /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
425: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);
426: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
427: MatCreateVecs(A,NULL,&b);
428: VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);
429: VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol);
430: ISDestroy(&is_iden);
431: VecDestroy(&b);
432: }
433: lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
434: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
436: PASTIX_CALL(&(lu->pastix_data),
437: lu->pastix_comm,
438: lu->n,
439: lu->colptr,
440: lu->row,
441: (PastixScalar*)lu->val,
442: lu->perm,
443: lu->invp,
444: (PastixScalar*)lu->rhs,
445: lu->rhsnbr,
446: lu->iparm,
447: lu->dparm);
448: 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]);
449: } else {
450: lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
451: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
452: PASTIX_CALL(&(lu->pastix_data),
453: lu->pastix_comm,
454: lu->n,
455: lu->colptr,
456: lu->row,
457: (PastixScalar*)lu->val,
458: lu->perm,
459: lu->invp,
460: (PastixScalar*)lu->rhs,
461: lu->rhsnbr,
462: lu->iparm,
463: lu->dparm);
464: 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]);
465: }
467: (F)->assembled = PETSC_TRUE;
468: lu->matstruc = SAME_NONZERO_PATTERN;
469: lu->CleanUpPastix = PETSC_TRUE;
470: return(0);
471: }
473: /* Note the Petsc r and c permutations are ignored */
474: PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
475: {
476: Mat_Pastix *lu = (Mat_Pastix*)F->data;
479: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
480: lu->iparm[IPARM_SYM] = API_SYM_YES;
481: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
482: F->ops->lufactornumeric = MatFactorNumeric_PaStiX;
483: return(0);
484: }
487: /* Note the Petsc r permutation is ignored */
488: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
489: {
490: Mat_Pastix *lu = (Mat_Pastix*)(F)->data;
493: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT;
494: lu->iparm[IPARM_SYM] = API_SYM_NO;
495: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
496: (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
497: return(0);
498: }
500: PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
501: {
502: PetscErrorCode ierr;
503: PetscBool iascii;
504: PetscViewerFormat format;
507: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
508: if (iascii) {
509: PetscViewerGetFormat(viewer,&format);
510: if (format == PETSC_VIEWER_ASCII_INFO) {
511: Mat_Pastix *lu=(Mat_Pastix*)A->data;
513: PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");
514: PetscViewerASCIIPrintf(viewer," Matrix type : %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
515: PetscViewerASCIIPrintf(viewer," Level of printing (0,1,2): %d \n",lu->iparm[IPARM_VERBOSE]);
516: PetscViewerASCIIPrintf(viewer," Number of refinements iterations : %d \n",lu->iparm[IPARM_NBITER]);
517: PetscPrintf(PETSC_COMM_SELF," Error : %g \n",lu->dparm[DPARM_RELATIVE_ERROR]);
518: }
519: }
520: return(0);
521: }
524: /*MC
525: MATSOLVERPASTIX - A solver package providing direct solvers (LU) for distributed
526: and sequential matrices via the external package PaStiX.
528: Use ./configure --download-pastix --download-ptscotch to have PETSc installed with PasTiX
530: Use -pc_type lu -pc_factor_mat_solver_package pastix to use this direct solver
532: Options Database Keys:
533: + -mat_pastix_verbose <0,1,2> - print level
534: - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task.
536: Level: beginner
538: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
540: M*/
543: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
544: {
545: Mat_Pastix *lu =(Mat_Pastix*)A->data;
548: info->block_size = 1.0;
549: info->nz_allocated = lu->iparm[IPARM_NNZEROS];
550: info->nz_used = lu->iparm[IPARM_NNZEROS];
551: info->nz_unneeded = 0.0;
552: info->assemblies = 0.0;
553: info->mallocs = 0.0;
554: info->memory = 0.0;
555: info->fill_ratio_given = 0;
556: info->fill_ratio_needed = 0;
557: info->factor_mallocs = 0;
558: return(0);
559: }
561: static PetscErrorCode MatFactorGetSolverPackage_pastix(Mat A,const MatSolverPackage *type)
562: {
564: *type = MATSOLVERPASTIX;
565: return(0);
566: }
568: /*
569: The seq and mpi versions of this function are the same
570: */
571: static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
572: {
573: Mat B;
575: Mat_Pastix *pastix;
578: if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
579: /* Create the factorization matrix */
580: MatCreate(PetscObjectComm((PetscObject)A),&B);
581: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
582: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
583: MatSetUp(B);
585: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
586: B->ops->view = MatView_PaStiX;
587: B->ops->getinfo = MatGetInfo_PaStiX;
589: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
591: B->factortype = MAT_FACTOR_LU;
592:
593: /* set solvertype */
594: PetscFree(B->solvertype);
595: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
597: PetscNewLog(B,&pastix);
599: pastix->CleanUpPastix = PETSC_FALSE;
600: pastix->scat_rhs = NULL;
601: pastix->scat_sol = NULL;
602: B->ops->getinfo = MatGetInfo_External;
603: B->ops->destroy = MatDestroy_Pastix;
604: B->data = (void*)pastix;
606: *F = B;
607: return(0);
608: }
610: static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
611: {
612: Mat B;
614: Mat_Pastix *pastix;
617: if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
618: /* Create the factorization matrix */
619: MatCreate(PetscObjectComm((PetscObject)A),&B);
620: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
621: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
622: MatSetUp(B);
624: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
625: B->ops->view = MatView_PaStiX;
626: B->ops->getinfo = MatGetInfo_PaStiX;
627: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
629: B->factortype = MAT_FACTOR_LU;
631: /* set solvertype */
632: PetscFree(B->solvertype);
633: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
635: PetscNewLog(B,&pastix);
637: pastix->CleanUpPastix = PETSC_FALSE;
638: pastix->scat_rhs = NULL;
639: pastix->scat_sol = NULL;
640: B->ops->getinfo = MatGetInfo_External;
641: B->ops->destroy = MatDestroy_Pastix;
642: B->data = (void*)pastix;
644: *F = B;
645: return(0);
646: }
648: static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
649: {
650: Mat B;
652: Mat_Pastix *pastix;
655: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
656: /* Create the factorization matrix */
657: MatCreate(PetscObjectComm((PetscObject)A),&B);
658: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
659: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
660: MatSetUp(B);
662: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
663: B->ops->view = MatView_PaStiX;
664: B->ops->getinfo = MatGetInfo_PaStiX;
665: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
667: B->factortype = MAT_FACTOR_CHOLESKY;
669: /* set solvertype */
670: PetscFree(B->solvertype);
671: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
673: PetscNewLog(B,&pastix);
675: pastix->CleanUpPastix = PETSC_FALSE;
676: pastix->scat_rhs = NULL;
677: pastix->scat_sol = NULL;
678: B->ops->getinfo = MatGetInfo_External;
679: B->ops->destroy = MatDestroy_Pastix;
680: B->data = (void*)pastix;
682: *F = B;
683: return(0);
684: }
686: static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
687: {
688: Mat B;
690: Mat_Pastix *pastix;
693: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
695: /* Create the factorization matrix */
696: MatCreate(PetscObjectComm((PetscObject)A),&B);
697: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
698: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
699: MatSetUp(B);
701: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
702: B->ops->view = MatView_PaStiX;
703: B->ops->getinfo = MatGetInfo_PaStiX;
704: B->ops->destroy = MatDestroy_Pastix;
705: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);
707: B->factortype = MAT_FACTOR_CHOLESKY;
709: /* set solvertype */
710: PetscFree(B->solvertype);
711: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
713: PetscNewLog(B,&pastix);
715: pastix->CleanUpPastix = PETSC_FALSE;
716: pastix->scat_rhs = NULL;
717: pastix->scat_sol = NULL;
718: B->data = (void*)pastix;
720: *F = B;
721: return(0);
722: }
724: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Pastix(void)
725: {
729: MatSolverPackageRegister(MATSOLVERPASTIX,MATMPIAIJ, MAT_FACTOR_LU,MatGetFactor_mpiaij_pastix);
730: MatSolverPackageRegister(MATSOLVERPASTIX,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_pastix);
731: MatSolverPackageRegister(MATSOLVERPASTIX,MATMPISBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_mpisbaij_pastix);
732: MatSolverPackageRegister(MATSOLVERPASTIX,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_pastix);
733: return(0);
734: }