Actual source code: pastix.c
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;
84: PetscInt colidx;
85: PetscInt *colcount;
86: PetscBool isSBAIJ;
87: PetscBool isSeqSBAIJ;
88: PetscBool isMpiSBAIJ;
89: PetscBool isSym;
91: MatIsSymmetric(A, 0.0, &isSym);
92: PetscObjectTypeCompare((PetscObject)A, MATSBAIJ, &isSBAIJ);
93: PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &isSeqSBAIJ);
94: PetscObjectTypeCompare((PetscObject)A, MATMPISBAIJ, &isMpiSBAIJ);
96: *n = A->cmap->N;
98: /* PaStiX only needs triangular matrix if matrix is symmetric
99: */
100: if (isSym && !(isSBAIJ || isSeqSBAIJ || isMpiSBAIJ)) nnz = (aa->nz - *n) / 2 + *n;
101: else nnz = aa->nz;
103: if (!valOnly) {
104: PetscMalloc1((*n) + 1, colptr);
105: PetscMalloc1(nnz, row);
106: PetscMalloc1(nnz, values);
108: if (isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) {
109: PetscArraycpy(*colptr, rowptr, (*n) + 1);
110: for (i = 0; i < *n + 1; i++) (*colptr)[i] += base;
111: PetscArraycpy(*row, col, nnz);
112: for (i = 0; i < nnz; i++) (*row)[i] += base;
113: PetscArraycpy(*values, rvalues, nnz);
114: } else {
115: PetscMalloc1(*n, &colcount);
117: for (i = 0; i < m; i++) colcount[i] = 0;
118: /* Fill-in colptr */
119: for (i = 0; i < m; i++) {
120: for (j = rowptr[i]; j < rowptr[i + 1]; j++) {
121: if (!isSym || col[j] <= i) colcount[col[j]]++;
122: }
123: }
125: (*colptr)[0] = base;
126: for (j = 0; j < *n; j++) {
127: (*colptr)[j + 1] = (*colptr)[j] + colcount[j];
128: /* in next loop we fill starting from (*colptr)[colidx] - base */
129: colcount[j] = -base;
130: }
132: /* Fill-in rows and values */
133: for (i = 0; i < m; i++) {
134: for (j = rowptr[i]; j < rowptr[i + 1]; j++) {
135: if (!isSym || col[j] <= i) {
136: colidx = col[j];
137: idx = (*colptr)[colidx] + colcount[colidx];
138: (*row)[idx] = i + base;
139: (*values)[idx] = rvalues[j];
140: colcount[colidx]++;
141: }
142: }
143: }
144: PetscFree(colcount);
145: }
146: } else {
147: /* Fill-in only values */
148: for (i = 0; i < m; i++) {
149: for (j = rowptr[i]; j < rowptr[i + 1]; j++) {
150: colidx = col[j];
151: if ((isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) || !isSym || col[j] <= i) {
152: /* look for the value to fill */
153: for (k = (*colptr)[colidx] - base; k < (*colptr)[colidx + 1] - base; k++) {
154: if (((*row)[k] - base) == i) {
155: (*values)[k] = rvalues[j];
156: break;
157: }
158: }
159: /* data structure of sparse matrix has changed */
161: }
162: }
163: }
164: }
165: return 0;
166: }
168: /*
169: Call clean step of PaStiX if lu->CleanUpPastix == true.
170: Free the CSC matrix.
171: */
172: PetscErrorCode MatDestroy_Pastix(Mat A)
173: {
174: Mat_Pastix *lu = (Mat_Pastix *)A->data;
176: if (lu->CleanUpPastix) {
177: /* Terminate instance, deallocate memories */
178: VecScatterDestroy(&lu->scat_rhs);
179: VecDestroy(&lu->b_seq);
180: VecScatterDestroy(&lu->scat_sol);
182: lu->iparm[IPARM_START_TASK] = API_TASK_CLEAN;
183: lu->iparm[IPARM_END_TASK] = API_TASK_CLEAN;
185: PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm);
187: PetscFree(lu->colptr);
188: PetscFree(lu->row);
189: PetscFree(lu->val);
190: PetscFree(lu->perm);
191: PetscFree(lu->invp);
192: MPI_Comm_free(&(lu->pastix_comm));
193: }
194: PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL);
195: PetscFree(A->data);
196: return 0;
197: }
199: /*
200: Gather right-hand-side.
201: Call for Solve step.
202: Scatter solution.
203: */
204: PetscErrorCode MatSolve_PaStiX(Mat A, Vec b, Vec x)
205: {
206: Mat_Pastix *lu = (Mat_Pastix *)A->data;
207: PetscScalar *array;
208: Vec x_seq;
210: lu->rhsnbr = 1;
211: x_seq = lu->b_seq;
212: if (lu->commSize > 1) {
213: /* PaStiX only supports centralized rhs. Scatter b into a sequential rhs vector */
214: VecScatterBegin(lu->scat_rhs, b, x_seq, INSERT_VALUES, SCATTER_FORWARD);
215: VecScatterEnd(lu->scat_rhs, b, x_seq, INSERT_VALUES, SCATTER_FORWARD);
216: VecGetArray(x_seq, &array);
217: } else { /* size == 1 */
218: VecCopy(b, x);
219: VecGetArray(x, &array);
220: }
221: lu->rhs = array;
222: if (lu->commSize == 1) {
223: VecRestoreArray(x, &array);
224: } else {
225: VecRestoreArray(x_seq, &array);
226: }
228: /* solve phase */
229: /*-------------*/
230: lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
231: lu->iparm[IPARM_END_TASK] = API_TASK_REFINE;
232: lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;
234: PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm);
237: if (lu->commSize == 1) {
238: VecRestoreArray(x, &(lu->rhs));
239: } else {
240: VecRestoreArray(x_seq, &(lu->rhs));
241: }
243: if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
244: VecScatterBegin(lu->scat_sol, x_seq, x, INSERT_VALUES, SCATTER_FORWARD);
245: VecScatterEnd(lu->scat_sol, x_seq, x, INSERT_VALUES, SCATTER_FORWARD);
246: }
247: return 0;
248: }
250: /*
251: Numeric factorisation using PaStiX solver.
253: */
254: PetscErrorCode MatFactorNumeric_PaStiX(Mat F, Mat A, const MatFactorInfo *info)
255: {
256: Mat_Pastix *lu = (Mat_Pastix *)(F)->data;
257: Mat *tseq;
258: PetscInt icntl;
259: PetscInt M = A->rmap->N;
260: PetscBool valOnly, flg, isSym;
261: IS is_iden;
262: Vec b;
263: IS isrow;
264: PetscBool isSeqAIJ, isSeqSBAIJ, isMPIAIJ;
266: PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSeqAIJ);
267: PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isMPIAIJ);
268: PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &isSeqSBAIJ);
269: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
270: (F)->ops->solve = MatSolve_PaStiX;
272: /* Initialize a PASTIX instance */
273: MPI_Comm_dup(PetscObjectComm((PetscObject)A), &(lu->pastix_comm));
274: MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
275: MPI_Comm_size(lu->pastix_comm, &lu->commSize);
277: /* Set pastix options */
278: lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
279: lu->iparm[IPARM_START_TASK] = API_TASK_INIT;
280: lu->iparm[IPARM_END_TASK] = API_TASK_INIT;
282: lu->rhsnbr = 1;
284: /* Call to set default pastix options */
285: PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm);
288: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "PaStiX Options", "Mat");
289: icntl = -1;
290: lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
291: PetscOptionsInt("-mat_pastix_verbose", "iparm[IPARM_VERBOSE] : level of printing (0 to 2)", "None", lu->iparm[IPARM_VERBOSE], &icntl, &flg);
292: if ((flg && icntl >= 0) || PetscLogPrintInfo) lu->iparm[IPARM_VERBOSE] = icntl;
293: icntl = -1;
294: PetscOptionsInt("-mat_pastix_threadnbr", "iparm[IPARM_THREAD_NBR] : Number of thread by MPI node", "None", lu->iparm[IPARM_THREAD_NBR], &icntl, &flg);
295: if ((flg && icntl > 0)) lu->iparm[IPARM_THREAD_NBR] = icntl;
296: PetscOptionsEnd();
297: valOnly = PETSC_FALSE;
298: } else {
299: if (isSeqAIJ || isMPIAIJ) {
300: PetscFree(lu->colptr);
301: PetscFree(lu->row);
302: PetscFree(lu->val);
303: valOnly = PETSC_FALSE;
304: } else valOnly = PETSC_TRUE;
305: }
307: lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES;
309: /* convert mpi A to seq mat A */
310: ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &isrow);
311: MatCreateSubMatrices(A, 1, &isrow, &isrow, MAT_INITIAL_MATRIX, &tseq);
312: ISDestroy(&isrow);
314: MatConvertToCSC(*tseq, valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
315: MatIsSymmetric(*tseq, 0.0, &isSym);
316: MatDestroyMatrices(1, &tseq);
318: if (!lu->perm) {
319: PetscMalloc1(lu->n, &(lu->perm));
320: PetscMalloc1(lu->n, &(lu->invp));
321: }
323: if (isSym) {
324: /* On symmetric matrix, LLT */
325: lu->iparm[IPARM_SYM] = API_SYM_YES;
326: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
327: } else {
328: /* On unsymmetric matrix, LU */
329: lu->iparm[IPARM_SYM] = API_SYM_NO;
330: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
331: }
333: /*----------------*/
334: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
335: if (!(isSeqAIJ || isSeqSBAIJ) && !lu->b_seq) {
336: /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
337: VecCreateSeq(PETSC_COMM_SELF, A->cmap->N, &lu->b_seq);
338: ISCreateStride(PETSC_COMM_SELF, A->cmap->N, 0, 1, &is_iden);
339: MatCreateVecs(A, NULL, &b);
340: VecScatterCreate(b, is_iden, lu->b_seq, is_iden, &lu->scat_rhs);
341: VecScatterCreate(lu->b_seq, is_iden, b, is_iden, &lu->scat_sol);
342: ISDestroy(&is_iden);
343: VecDestroy(&b);
344: }
345: lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
346: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
348: PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm);
350: } else {
351: lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
352: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
353: PASTIX_CALL(&(lu->pastix_data), lu->pastix_comm, lu->n, lu->colptr, lu->row, (PastixScalar *)lu->val, lu->perm, lu->invp, (PastixScalar *)lu->rhs, lu->rhsnbr, lu->iparm, lu->dparm);
355: }
357: (F)->assembled = PETSC_TRUE;
358: lu->matstruc = SAME_NONZERO_PATTERN;
359: lu->CleanUpPastix = PETSC_TRUE;
360: return 0;
361: }
363: /* Note the Petsc r and c permutations are ignored */
364: PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
365: {
366: Mat_Pastix *lu = (Mat_Pastix *)F->data;
368: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
369: lu->iparm[IPARM_SYM] = API_SYM_YES;
370: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
371: F->ops->lufactornumeric = MatFactorNumeric_PaStiX;
372: return 0;
373: }
375: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F, Mat A, IS r, const MatFactorInfo *info)
376: {
377: Mat_Pastix *lu = (Mat_Pastix *)(F)->data;
379: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT;
380: lu->iparm[IPARM_SYM] = API_SYM_NO;
381: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
382: (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
383: return 0;
384: }
386: PetscErrorCode MatView_PaStiX(Mat A, PetscViewer viewer)
387: {
388: PetscBool iascii;
389: PetscViewerFormat format;
391: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
392: if (iascii) {
393: PetscViewerGetFormat(viewer, &format);
394: if (format == PETSC_VIEWER_ASCII_INFO) {
395: Mat_Pastix *lu = (Mat_Pastix *)A->data;
397: PetscViewerASCIIPrintf(viewer, "PaStiX run parameters:\n");
398: PetscViewerASCIIPrintf(viewer, " Matrix type : %s \n", ((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
399: PetscViewerASCIIPrintf(viewer, " Level of printing (0,1,2): %" PetscInt_FMT " \n", lu->iparm[IPARM_VERBOSE]);
400: PetscViewerASCIIPrintf(viewer, " Number of refinements iterations : %" PetscInt_FMT " \n", lu->iparm[IPARM_NBITER]);
401: PetscPrintf(PETSC_COMM_SELF, " Error : %g \n", lu->dparm[DPARM_RELATIVE_ERROR]);
402: }
403: }
404: return 0;
405: }
407: /*MC
408: MATSOLVERPASTIX - A solver package providing direct solvers (LU) for distributed
409: and sequential matrices via the external package PaStiX.
411: Use ./configure --download-pastix --download-ptscotch to have PETSc installed with PasTiX
413: Use -pc_type lu -pc_factor_mat_solver_type pastix to use this direct solver
415: Options Database Keys:
416: + -mat_pastix_verbose <0,1,2> - print level
417: - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task.
419: Notes:
420: This only works for matrices with symmetric nonzero structure, if you pass it a matrix with
421: nonsymmetric structure PasTiX and hence PETSc return with an error.
423: Level: beginner
425: .seealso: `PCFactorSetMatSolverType()`, `MatSolverType`, `MatGetFactor()`
426: M*/
428: PetscErrorCode MatGetInfo_PaStiX(Mat A, MatInfoType flag, MatInfo *info)
429: {
430: Mat_Pastix *lu = (Mat_Pastix *)A->data;
432: info->block_size = 1.0;
433: info->nz_allocated = lu->iparm[IPARM_NNZEROS];
434: info->nz_used = lu->iparm[IPARM_NNZEROS];
435: info->nz_unneeded = 0.0;
436: info->assemblies = 0.0;
437: info->mallocs = 0.0;
438: info->memory = 0.0;
439: info->fill_ratio_given = 0;
440: info->fill_ratio_needed = 0;
441: info->factor_mallocs = 0;
442: return 0;
443: }
445: static PetscErrorCode MatFactorGetSolverType_pastix(Mat A, MatSolverType *type)
446: {
447: *type = MATSOLVERPASTIX;
448: return 0;
449: }
451: /*
452: The seq and mpi versions of this function are the same
453: */
454: static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A, MatFactorType ftype, Mat *F)
455: {
456: Mat B;
457: Mat_Pastix *pastix;
460: /* Create the factorization matrix */
461: MatCreate(PetscObjectComm((PetscObject)A), &B);
462: MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
463: PetscStrallocpy("pastix", &((PetscObject)B)->type_name);
464: MatSetUp(B);
466: B->trivialsymbolic = PETSC_TRUE;
467: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
468: B->ops->view = MatView_PaStiX;
469: B->ops->getinfo = MatGetInfo_PaStiX;
471: PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_pastix);
473: B->factortype = MAT_FACTOR_LU;
475: /* set solvertype */
476: PetscFree(B->solvertype);
477: PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype);
479: PetscNew(&pastix);
481: pastix->CleanUpPastix = PETSC_FALSE;
482: pastix->scat_rhs = NULL;
483: pastix->scat_sol = NULL;
484: B->ops->getinfo = MatGetInfo_External;
485: B->ops->destroy = MatDestroy_Pastix;
486: B->data = (void *)pastix;
488: *F = B;
489: return 0;
490: }
492: static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A, MatFactorType ftype, Mat *F)
493: {
494: Mat B;
495: Mat_Pastix *pastix;
498: /* Create the factorization matrix */
499: MatCreate(PetscObjectComm((PetscObject)A), &B);
500: MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
501: PetscStrallocpy("pastix", &((PetscObject)B)->type_name);
502: MatSetUp(B);
504: B->trivialsymbolic = PETSC_TRUE;
505: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
506: B->ops->view = MatView_PaStiX;
507: B->ops->getinfo = MatGetInfo_PaStiX;
508: PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_pastix);
510: B->factortype = MAT_FACTOR_LU;
512: /* set solvertype */
513: PetscFree(B->solvertype);
514: PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype);
516: PetscNew(&pastix);
518: pastix->CleanUpPastix = PETSC_FALSE;
519: pastix->scat_rhs = NULL;
520: pastix->scat_sol = NULL;
521: B->ops->getinfo = MatGetInfo_External;
522: B->ops->destroy = MatDestroy_Pastix;
523: B->data = (void *)pastix;
525: *F = B;
526: return 0;
527: }
529: static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A, MatFactorType ftype, Mat *F)
530: {
531: Mat B;
532: Mat_Pastix *pastix;
535: /* Create the factorization matrix */
536: MatCreate(PetscObjectComm((PetscObject)A), &B);
537: MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
538: PetscStrallocpy("pastix", &((PetscObject)B)->type_name);
539: MatSetUp(B);
541: B->trivialsymbolic = PETSC_TRUE;
542: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
543: B->ops->view = MatView_PaStiX;
544: B->ops->getinfo = MatGetInfo_PaStiX;
545: PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_pastix);
547: B->factortype = MAT_FACTOR_CHOLESKY;
549: /* set solvertype */
550: PetscFree(B->solvertype);
551: PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype);
553: PetscNew(&pastix);
555: pastix->CleanUpPastix = PETSC_FALSE;
556: pastix->scat_rhs = NULL;
557: pastix->scat_sol = NULL;
558: B->ops->getinfo = MatGetInfo_External;
559: B->ops->destroy = MatDestroy_Pastix;
560: B->data = (void *)pastix;
561: *F = B;
562: return 0;
563: }
565: static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A, MatFactorType ftype, Mat *F)
566: {
567: Mat B;
568: Mat_Pastix *pastix;
572: /* Create the factorization matrix */
573: MatCreate(PetscObjectComm((PetscObject)A), &B);
574: MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
575: PetscStrallocpy("pastix", &((PetscObject)B)->type_name);
576: MatSetUp(B);
578: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
579: B->ops->view = MatView_PaStiX;
580: B->ops->getinfo = MatGetInfo_PaStiX;
581: B->ops->destroy = MatDestroy_Pastix;
582: PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_pastix);
584: B->factortype = MAT_FACTOR_CHOLESKY;
586: /* set solvertype */
587: PetscFree(B->solvertype);
588: PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype);
590: PetscNew(&pastix);
592: pastix->CleanUpPastix = PETSC_FALSE;
593: pastix->scat_rhs = NULL;
594: pastix->scat_sol = NULL;
595: B->data = (void *)pastix;
597: *F = B;
598: return 0;
599: }
601: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Pastix(void)
602: {
603: MatSolverTypeRegister(MATSOLVERPASTIX, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_pastix);
604: MatSolverTypeRegister(MATSOLVERPASTIX, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_pastix);
605: MatSolverTypeRegister(MATSOLVERPASTIX, MATMPISBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_mpisbaij_pastix);
606: MatSolverTypeRegister(MATSOLVERPASTIX, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_pastix);
607: return 0;
608: }