Actual source code: pastix.c
petsc-3.11.4 2019-09-28
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: #else
21: #define PASTIX_CALL z_pastix
22: #endif
24: #else /* PETSC_USE_COMPLEX */
26: #if defined(PETSC_USE_REAL_SINGLE)
27: #define PASTIX_CALL s_pastix
28: #else
29: #define PASTIX_CALL d_pastix
30: #endif
32: #endif /* PETSC_USE_COMPLEX */
34: typedef PetscScalar PastixScalar;
36: typedef struct Mat_Pastix_ {
37: pastix_data_t *pastix_data; /* Pastix data storage structure */
38: MatStructure matstruc;
39: PetscInt n; /* Number of columns in the matrix */
40: PetscInt *colptr; /* Index of first element of each column in row and val */
41: PetscInt *row; /* Row of each element of the matrix */
42: PetscScalar *val; /* Value of each element of the matrix */
43: PetscInt *perm; /* Permutation tabular */
44: PetscInt *invp; /* Reverse permutation tabular */
45: PetscScalar *rhs; /* Rhight-hand-side member */
46: PetscInt rhsnbr; /* Rhight-hand-side number (must be 1) */
47: PetscInt iparm[IPARM_SIZE]; /* Integer parameters */
48: double dparm[DPARM_SIZE]; /* Floating point parameters */
49: MPI_Comm pastix_comm; /* PaStiX MPI communicator */
50: PetscMPIInt commRank; /* MPI rank */
51: PetscMPIInt commSize; /* MPI communicator size */
52: PetscBool CleanUpPastix; /* Boolean indicating if we call PaStiX clean step */
53: VecScatter scat_rhs;
54: VecScatter scat_sol;
55: Vec b_seq;
56: } Mat_Pastix;
58: extern PetscErrorCode MatDuplicate_Pastix(Mat,MatDuplicateOption,Mat*);
60: /*
61: convert Petsc seqaij matrix to CSC: colptr[n], row[nz], val[nz]
63: input:
64: A - matrix in seqaij or mpisbaij (bs=1) format
65: valOnly - FALSE: spaces are allocated and values are set for the CSC
66: TRUE: Only fill values
67: output:
68: n - Size of the matrix
69: colptr - Index of first element of each column in row and val
70: row - Row of each element of the matrix
71: values - Value of each element of the matrix
72: */
73: PetscErrorCode MatConvertToCSC(Mat A,PetscBool valOnly,PetscInt *n,PetscInt **colptr,PetscInt **row,PetscScalar **values)
74: {
75: Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data;
76: PetscInt *rowptr = aa->i;
77: PetscInt *col = aa->j;
78: PetscScalar *rvalues = aa->a;
79: PetscInt m = A->rmap->N;
80: PetscInt nnz;
81: PetscInt i,j, k;
82: PetscInt base = 1;
83: PetscInt idx;
85: PetscInt colidx;
86: PetscInt *colcount;
87: PetscBool isSBAIJ;
88: PetscBool isSeqSBAIJ;
89: PetscBool isMpiSBAIJ;
90: PetscBool isSym;
93: MatIsSymmetric(A,0.0,&isSym);
94: PetscObjectTypeCompare((PetscObject)A,MATSBAIJ,&isSBAIJ);
95: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
96: PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMpiSBAIJ);
98: *n = A->cmap->N;
100: /* PaStiX only needs triangular matrix if matrix is symmetric
101: */
102: if (isSym && !(isSBAIJ || isSeqSBAIJ || isMpiSBAIJ)) nnz = (aa->nz - *n)/2 + *n;
103: else nnz = aa->nz;
105: if (!valOnly) {
106: PetscMalloc1((*n)+1,colptr);
107: PetscMalloc1(nnz,row);
108: PetscMalloc1(nnz,values);
110: if (isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) {
111: PetscMemcpy (*colptr, rowptr, ((*n)+1)*sizeof(PetscInt));
112: for (i = 0; i < *n+1; i++) (*colptr)[i] += base;
113: PetscMemcpy (*row, col, (nnz)*sizeof(PetscInt));
114: for (i = 0; i < nnz; i++) (*row)[i] += base;
115: PetscMemcpy (*values, rvalues, (nnz)*sizeof(PetscScalar));
116: } else {
117: PetscMalloc1(*n,&colcount);
119: for (i = 0; i < m; i++) colcount[i] = 0;
120: /* Fill-in colptr */
121: for (i = 0; i < m; i++) {
122: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
123: if (!isSym || col[j] <= i) colcount[col[j]]++;
124: }
125: }
127: (*colptr)[0] = base;
128: for (j = 0; j < *n; j++) {
129: (*colptr)[j+1] = (*colptr)[j] + colcount[j];
130: /* in next loop we fill starting from (*colptr)[colidx] - base */
131: colcount[j] = -base;
132: }
134: /* Fill-in rows and values */
135: for (i = 0; i < m; i++) {
136: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
137: if (!isSym || col[j] <= i) {
138: colidx = col[j];
139: idx = (*colptr)[colidx] + colcount[colidx];
140: (*row)[idx] = i + base;
141: (*values)[idx] = rvalues[j];
142: colcount[colidx]++;
143: }
144: }
145: }
146: PetscFree(colcount);
147: }
148: } else {
149: /* Fill-in only values */
150: for (i = 0; i < m; i++) {
151: for (j = rowptr[i]; j < rowptr[i+1]; j++) {
152: colidx = col[j];
153: if ((isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) ||!isSym || col[j] <= i) {
154: /* look for the value to fill */
155: for (k = (*colptr)[colidx] - base; k < (*colptr)[colidx + 1] - base; k++) {
156: if (((*row)[k]-base) == i) {
157: (*values)[k] = rvalues[j];
158: break;
159: }
160: }
161: /* data structure of sparse matrix has changed */
162: if (k == (*colptr)[colidx + 1] - base) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"overflow on k %D",k);
163: }
164: }
165: }
166: }
167: return(0);
168: }
170: /*
171: Call clean step of PaStiX if lu->CleanUpPastix == true.
172: Free the CSC matrix.
173: */
174: PetscErrorCode MatDestroy_Pastix(Mat A)
175: {
176: Mat_Pastix *lu=(Mat_Pastix*)A->data;
180: if (lu->CleanUpPastix) {
181: /* Terminate instance, deallocate memories */
182: VecScatterDestroy(&lu->scat_rhs);
183: VecDestroy(&lu->b_seq);
184: VecScatterDestroy(&lu->scat_sol);
186: lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN;
187: lu->iparm[IPARM_END_TASK] =API_TASK_CLEAN;
189: PASTIX_CALL(&(lu->pastix_data),
190: lu->pastix_comm,
191: lu->n,
192: lu->colptr,
193: lu->row,
194: (PastixScalar*)lu->val,
195: lu->perm,
196: lu->invp,
197: (PastixScalar*)lu->rhs,
198: lu->rhsnbr,
199: lu->iparm,
200: lu->dparm);
201: if (lu->iparm[IPARM_ERROR_NUMBER] != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in destroy: iparm(IPARM_ERROR_NUMBER)=%d\n",lu->iparm[IPARM_ERROR_NUMBER]);
202: PetscFree(lu->colptr);
203: PetscFree(lu->row);
204: PetscFree(lu->val);
205: PetscFree(lu->perm);
206: PetscFree(lu->invp);
207: MPI_Comm_free(&(lu->pastix_comm));
208: }
209: PetscFree(A->data);
210: return(0);
211: }
213: /*
214: Gather right-hand-side.
215: Call for Solve step.
216: Scatter solution.
217: */
218: PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
219: {
220: Mat_Pastix *lu=(Mat_Pastix*)A->data;
221: PetscScalar *array;
222: Vec x_seq;
226: lu->rhsnbr = 1;
227: x_seq = lu->b_seq;
228: if (lu->commSize > 1) {
229: /* PaStiX only supports centralized rhs. Scatter b into a seqential rhs vector */
230: VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
231: VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
232: VecGetArray(x_seq,&array);
233: } else { /* size == 1 */
234: VecCopy(b,x);
235: VecGetArray(x,&array);
236: }
237: lu->rhs = array;
238: if (lu->commSize == 1) {
239: VecRestoreArray(x,&array);
240: } else {
241: VecRestoreArray(x_seq,&array);
242: }
244: /* solve phase */
245: /*-------------*/
246: lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
247: lu->iparm[IPARM_END_TASK] = API_TASK_REFINE;
248: lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;
250: PASTIX_CALL(&(lu->pastix_data),
251: lu->pastix_comm,
252: lu->n,
253: lu->colptr,
254: lu->row,
255: (PastixScalar*)lu->val,
256: lu->perm,
257: lu->invp,
258: (PastixScalar*)lu->rhs,
259: lu->rhsnbr,
260: lu->iparm,
261: lu->dparm);
262: 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]);
264: if (lu->commSize == 1) {
265: VecRestoreArray(x,&(lu->rhs));
266: } else {
267: VecRestoreArray(x_seq,&(lu->rhs));
268: }
270: if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
271: VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
272: VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
273: }
274: return(0);
275: }
277: /*
278: Numeric factorisation using PaStiX solver.
280: */
281: PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
282: {
283: Mat_Pastix *lu =(Mat_Pastix*)(F)->data;
284: Mat *tseq;
285: PetscErrorCode 0;
286: PetscInt icntl;
287: PetscInt M=A->rmap->N;
288: PetscBool valOnly,flg, isSym;
289: IS is_iden;
290: Vec b;
291: IS isrow;
292: PetscBool isSeqAIJ,isSeqSBAIJ,isMPIAIJ;
295: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
296: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
297: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
298: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
299: (F)->ops->solve = MatSolve_PaStiX;
301: /* Initialize a PASTIX instance */
302: MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm));
303: MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
304: MPI_Comm_size(lu->pastix_comm, &lu->commSize);
306: /* Set pastix options */
307: lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
308: lu->iparm[IPARM_START_TASK] = API_TASK_INIT;
309: lu->iparm[IPARM_END_TASK] = API_TASK_INIT;
311: lu->rhsnbr = 1;
313: /* Call to set default pastix options */
314: PASTIX_CALL(&(lu->pastix_data),
315: lu->pastix_comm,
316: lu->n,
317: lu->colptr,
318: lu->row,
319: (PastixScalar*)lu->val,
320: lu->perm,
321: lu->invp,
322: (PastixScalar*)lu->rhs,
323: lu->rhsnbr,
324: lu->iparm,
325: lu->dparm);
326: if (lu->iparm[IPARM_ERROR_NUMBER] != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by PaStiX in MatFactorNumeric: iparm(IPARM_ERROR_NUMBER)=%d\n",lu->iparm[IPARM_ERROR_NUMBER]);
328: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"PaStiX Options","Mat");
329: icntl = -1;
330: lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
331: PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg);
332: if ((flg && icntl >= 0) || PetscLogPrintInfo) {
333: lu->iparm[IPARM_VERBOSE] = icntl;
334: }
335: icntl=-1;
336: PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg);
337: if ((flg && icntl > 0)) {
338: lu->iparm[IPARM_THREAD_NBR] = icntl;
339: }
340: PetscOptionsEnd();
341: valOnly = PETSC_FALSE;
342: } else {
343: if (isSeqAIJ || isMPIAIJ) {
344: PetscFree(lu->colptr);
345: PetscFree(lu->row);
346: PetscFree(lu->val);
347: valOnly = PETSC_FALSE;
348: } else valOnly = PETSC_TRUE;
349: }
351: lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES;
353: /* convert mpi A to seq mat A */
354: ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);
355: MatCreateSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);
356: ISDestroy(&isrow);
358: MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
359: MatIsSymmetric(*tseq,0.0,&isSym);
360: MatDestroyMatrices(1,&tseq);
362: if (!lu->perm) {
363: PetscMalloc1(lu->n,&(lu->perm));
364: PetscMalloc1(lu->n,&(lu->invp));
365: }
367: if (isSym) {
368: /* On symmetric matrix, LLT */
369: lu->iparm[IPARM_SYM] = API_SYM_YES;
370: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
371: } else {
372: /* On unsymmetric matrix, LU */
373: lu->iparm[IPARM_SYM] = API_SYM_NO;
374: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
375: }
377: /*----------------*/
378: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
379: if (!(isSeqAIJ || isSeqSBAIJ) && !lu->b_seq) {
380: /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
381: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);
382: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
383: MatCreateVecs(A,NULL,&b);
384: VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);
385: VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol);
386: ISDestroy(&is_iden);
387: VecDestroy(&b);
388: }
389: lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
390: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
392: PASTIX_CALL(&(lu->pastix_data),
393: lu->pastix_comm,
394: lu->n,
395: lu->colptr,
396: lu->row,
397: (PastixScalar*)lu->val,
398: lu->perm,
399: lu->invp,
400: (PastixScalar*)lu->rhs,
401: lu->rhsnbr,
402: lu->iparm,
403: lu->dparm);
404: 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]);
405: } else {
406: lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
407: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
408: PASTIX_CALL(&(lu->pastix_data),
409: lu->pastix_comm,
410: lu->n,
411: lu->colptr,
412: lu->row,
413: (PastixScalar*)lu->val,
414: lu->perm,
415: lu->invp,
416: (PastixScalar*)lu->rhs,
417: lu->rhsnbr,
418: lu->iparm,
419: lu->dparm);
420: 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]);
421: }
423: (F)->assembled = PETSC_TRUE;
424: lu->matstruc = SAME_NONZERO_PATTERN;
425: lu->CleanUpPastix = PETSC_TRUE;
426: return(0);
427: }
429: /* Note the Petsc r and c permutations are ignored */
430: PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
431: {
432: Mat_Pastix *lu = (Mat_Pastix*)F->data;
435: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
436: lu->iparm[IPARM_SYM] = API_SYM_YES;
437: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
438: F->ops->lufactornumeric = MatFactorNumeric_PaStiX;
439: return(0);
440: }
443: /* Note the Petsc r permutation is ignored */
444: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
445: {
446: Mat_Pastix *lu = (Mat_Pastix*)(F)->data;
449: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT;
450: lu->iparm[IPARM_SYM] = API_SYM_NO;
451: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
452: (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
453: return(0);
454: }
456: PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
457: {
458: PetscErrorCode ierr;
459: PetscBool iascii;
460: PetscViewerFormat format;
463: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
464: if (iascii) {
465: PetscViewerGetFormat(viewer,&format);
466: if (format == PETSC_VIEWER_ASCII_INFO) {
467: Mat_Pastix *lu=(Mat_Pastix*)A->data;
469: PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");
470: PetscViewerASCIIPrintf(viewer," Matrix type : %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
471: PetscViewerASCIIPrintf(viewer," Level of printing (0,1,2): %d \n",lu->iparm[IPARM_VERBOSE]);
472: PetscViewerASCIIPrintf(viewer," Number of refinements iterations : %d \n",lu->iparm[IPARM_NBITER]);
473: PetscPrintf(PETSC_COMM_SELF," Error : %g \n",lu->dparm[DPARM_RELATIVE_ERROR]);
474: }
475: }
476: return(0);
477: }
480: /*MC
481: MATSOLVERPASTIX - A solver package providing direct solvers (LU) for distributed
482: and sequential matrices via the external package PaStiX.
484: Use ./configure --download-pastix --download-ptscotch to have PETSc installed with PasTiX
486: Use -pc_type lu -pc_factor_mat_solver_type pastix to use this direct solver
488: Options Database Keys:
489: + -mat_pastix_verbose <0,1,2> - print level
490: - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task.
492: Notes:
493: This only works for matrices with symmetric nonzero structure, if you pass it a matrix with
494: nonsymmetric structure PasTiX and hence PETSc return with an error.
496: Level: beginner
498: .seealso: PCFactorSetMatSolverType(), MatSolverType
500: M*/
503: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
504: {
505: Mat_Pastix *lu =(Mat_Pastix*)A->data;
508: info->block_size = 1.0;
509: info->nz_allocated = lu->iparm[IPARM_NNZEROS];
510: info->nz_used = lu->iparm[IPARM_NNZEROS];
511: info->nz_unneeded = 0.0;
512: info->assemblies = 0.0;
513: info->mallocs = 0.0;
514: info->memory = 0.0;
515: info->fill_ratio_given = 0;
516: info->fill_ratio_needed = 0;
517: info->factor_mallocs = 0;
518: return(0);
519: }
521: static PetscErrorCode MatFactorGetSolverType_pastix(Mat A,MatSolverType *type)
522: {
524: *type = MATSOLVERPASTIX;
525: return(0);
526: }
528: /*
529: The seq and mpi versions of this function are the same
530: */
531: static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
532: {
533: Mat B;
535: Mat_Pastix *pastix;
538: if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
539: /* Create the factorization matrix */
540: MatCreate(PetscObjectComm((PetscObject)A),&B);
541: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
542: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
543: MatSetUp(B);
545: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
546: B->ops->view = MatView_PaStiX;
547: B->ops->getinfo = MatGetInfo_PaStiX;
549: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);
551: B->factortype = MAT_FACTOR_LU;
552:
553: /* set solvertype */
554: PetscFree(B->solvertype);
555: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
557: PetscNewLog(B,&pastix);
559: pastix->CleanUpPastix = PETSC_FALSE;
560: pastix->scat_rhs = NULL;
561: pastix->scat_sol = NULL;
562: B->ops->getinfo = MatGetInfo_External;
563: B->ops->destroy = MatDestroy_Pastix;
564: B->data = (void*)pastix;
566: *F = B;
567: return(0);
568: }
570: static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
571: {
572: Mat B;
574: Mat_Pastix *pastix;
577: if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
578: /* Create the factorization matrix */
579: MatCreate(PetscObjectComm((PetscObject)A),&B);
580: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
581: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
582: MatSetUp(B);
584: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
585: B->ops->view = MatView_PaStiX;
586: B->ops->getinfo = MatGetInfo_PaStiX;
587: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);
589: B->factortype = MAT_FACTOR_LU;
591: /* set solvertype */
592: PetscFree(B->solvertype);
593: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
595: PetscNewLog(B,&pastix);
597: pastix->CleanUpPastix = PETSC_FALSE;
598: pastix->scat_rhs = NULL;
599: pastix->scat_sol = NULL;
600: B->ops->getinfo = MatGetInfo_External;
601: B->ops->destroy = MatDestroy_Pastix;
602: B->data = (void*)pastix;
604: *F = B;
605: return(0);
606: }
608: static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
609: {
610: Mat B;
612: Mat_Pastix *pastix;
615: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
616: /* Create the factorization matrix */
617: MatCreate(PetscObjectComm((PetscObject)A),&B);
618: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
619: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
620: MatSetUp(B);
622: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
623: B->ops->view = MatView_PaStiX;
624: B->ops->getinfo = MatGetInfo_PaStiX;
625: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);
627: B->factortype = MAT_FACTOR_CHOLESKY;
629: /* set solvertype */
630: PetscFree(B->solvertype);
631: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
633: PetscNewLog(B,&pastix);
635: pastix->CleanUpPastix = PETSC_FALSE;
636: pastix->scat_rhs = NULL;
637: pastix->scat_sol = NULL;
638: B->ops->getinfo = MatGetInfo_External;
639: B->ops->destroy = MatDestroy_Pastix;
640: B->data = (void*)pastix;
642: *F = B;
643: return(0);
644: }
646: static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
647: {
648: Mat B;
650: Mat_Pastix *pastix;
653: if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
655: /* Create the factorization matrix */
656: MatCreate(PetscObjectComm((PetscObject)A),&B);
657: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
658: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
659: MatSetUp(B);
661: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
662: B->ops->view = MatView_PaStiX;
663: B->ops->getinfo = MatGetInfo_PaStiX;
664: B->ops->destroy = MatDestroy_Pastix;
665: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_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->data = (void*)pastix;
680: *F = B;
681: return(0);
682: }
684: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Pastix(void)
685: {
689: MatSolverTypeRegister(MATSOLVERPASTIX,MATMPIAIJ, MAT_FACTOR_LU,MatGetFactor_mpiaij_pastix);
690: MatSolverTypeRegister(MATSOLVERPASTIX,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_pastix);
691: MatSolverTypeRegister(MATSOLVERPASTIX,MATMPISBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_mpisbaij_pastix);
692: MatSolverTypeRegister(MATSOLVERPASTIX,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_pastix);
693: return(0);
694: }