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),
186: lu->pastix_comm,
187: lu->n,
188: lu->colptr,
189: lu->row,
190: (PastixScalar*)lu->val,
191: lu->perm,
192: lu->invp,
193: (PastixScalar*)lu->rhs,
194: lu->rhsnbr,
195: lu->iparm,
196: lu->dparm);
198: PetscFree(lu->colptr);
199: PetscFree(lu->row);
200: PetscFree(lu->val);
201: PetscFree(lu->perm);
202: PetscFree(lu->invp);
203: MPI_Comm_free(&(lu->pastix_comm));
204: }
205: PetscFree(A->data);
206: return 0;
207: }
209: /*
210: Gather right-hand-side.
211: Call for Solve step.
212: Scatter solution.
213: */
214: PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
215: {
216: Mat_Pastix *lu = (Mat_Pastix*)A->data;
217: PetscScalar *array;
218: Vec x_seq;
220: lu->rhsnbr = 1;
221: x_seq = lu->b_seq;
222: if (lu->commSize > 1) {
223: /* PaStiX only supports centralized rhs. Scatter b into a sequential rhs vector */
224: VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
225: VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
226: VecGetArray(x_seq,&array);
227: } else { /* size == 1 */
228: VecCopy(b,x);
229: VecGetArray(x,&array);
230: }
231: lu->rhs = array;
232: if (lu->commSize == 1) {
233: VecRestoreArray(x,&array);
234: } else {
235: VecRestoreArray(x_seq,&array);
236: }
238: /* solve phase */
239: /*-------------*/
240: lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
241: lu->iparm[IPARM_END_TASK] = API_TASK_REFINE;
242: lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;
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);
258: if (lu->commSize == 1) {
259: VecRestoreArray(x,&(lu->rhs));
260: } else {
261: VecRestoreArray(x_seq,&(lu->rhs));
262: }
264: if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
265: VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
266: VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
267: }
268: return 0;
269: }
271: /*
272: Numeric factorisation using PaStiX solver.
274: */
275: PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
276: {
277: Mat_Pastix *lu =(Mat_Pastix*)(F)->data;
278: Mat *tseq;
280: PetscInt icntl;
281: PetscInt M=A->rmap->N;
282: PetscBool valOnly,flg, isSym;
283: IS is_iden;
284: Vec b;
285: IS isrow;
286: PetscBool isSeqAIJ,isSeqSBAIJ,isMPIAIJ;
288: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
289: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
290: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
291: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
292: (F)->ops->solve = MatSolve_PaStiX;
294: /* Initialize a PASTIX instance */
295: MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm));
296: MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
297: MPI_Comm_size(lu->pastix_comm, &lu->commSize);
299: /* Set pastix options */
300: lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
301: lu->iparm[IPARM_START_TASK] = API_TASK_INIT;
302: lu->iparm[IPARM_END_TASK] = API_TASK_INIT;
304: lu->rhsnbr = 1;
306: /* Call to set default pastix options */
307: PASTIX_CALL(&(lu->pastix_data),
308: lu->pastix_comm,
309: lu->n,
310: lu->colptr,
311: lu->row,
312: (PastixScalar*)lu->val,
313: lu->perm,
314: lu->invp,
315: (PastixScalar*)lu->rhs,
316: lu->rhsnbr,
317: lu->iparm,
318: lu->dparm);
321: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"PaStiX Options","Mat");
322: icntl = -1;
323: lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;
324: PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg);
325: if ((flg && icntl >= 0) || PetscLogPrintInfo) {
326: lu->iparm[IPARM_VERBOSE] = icntl;
327: }
328: icntl=-1;
329: PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg);
330: if ((flg && icntl > 0)) {
331: lu->iparm[IPARM_THREAD_NBR] = icntl;
332: }
333: PetscOptionsEnd();
334: valOnly = PETSC_FALSE;
335: } else {
336: if (isSeqAIJ || isMPIAIJ) {
337: PetscFree(lu->colptr);
338: PetscFree(lu->row);
339: PetscFree(lu->val);
340: valOnly = PETSC_FALSE;
341: } else valOnly = PETSC_TRUE;
342: }
344: lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES;
346: /* convert mpi A to seq mat A */
347: ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);
348: MatCreateSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);
349: ISDestroy(&isrow);
351: MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
352: MatIsSymmetric(*tseq,0.0,&isSym);
353: MatDestroyMatrices(1,&tseq);
355: if (!lu->perm) {
356: PetscMalloc1(lu->n,&(lu->perm));
357: PetscMalloc1(lu->n,&(lu->invp));
358: }
360: if (isSym) {
361: /* On symmetric matrix, LLT */
362: lu->iparm[IPARM_SYM] = API_SYM_YES;
363: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
364: } else {
365: /* On unsymmetric matrix, LU */
366: lu->iparm[IPARM_SYM] = API_SYM_NO;
367: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
368: }
370: /*----------------*/
371: if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
372: if (!(isSeqAIJ || isSeqSBAIJ) && !lu->b_seq) {
373: /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
374: VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);
375: ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
376: MatCreateVecs(A,NULL,&b);
377: VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);
378: VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol);
379: ISDestroy(&is_iden);
380: VecDestroy(&b);
381: }
382: lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
383: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
385: PASTIX_CALL(&(lu->pastix_data),
386: lu->pastix_comm,
387: lu->n,
388: lu->colptr,
389: lu->row,
390: (PastixScalar*)lu->val,
391: lu->perm,
392: lu->invp,
393: (PastixScalar*)lu->rhs,
394: lu->rhsnbr,
395: lu->iparm,
396: lu->dparm);
398: } else {
399: lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
400: lu->iparm[IPARM_END_TASK] = API_TASK_NUMFACT;
401: PASTIX_CALL(&(lu->pastix_data),
402: lu->pastix_comm,
403: lu->n,
404: lu->colptr,
405: lu->row,
406: (PastixScalar*)lu->val,
407: lu->perm,
408: lu->invp,
409: (PastixScalar*)lu->rhs,
410: lu->rhsnbr,
411: lu->iparm,
412: lu->dparm);
414: }
416: (F)->assembled = PETSC_TRUE;
417: lu->matstruc = SAME_NONZERO_PATTERN;
418: lu->CleanUpPastix = PETSC_TRUE;
419: return 0;
420: }
422: /* Note the Petsc r and c permutations are ignored */
423: PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
424: {
425: Mat_Pastix *lu = (Mat_Pastix*)F->data;
427: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
428: lu->iparm[IPARM_SYM] = API_SYM_YES;
429: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
430: F->ops->lufactornumeric = MatFactorNumeric_PaStiX;
431: return 0;
432: }
434: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
435: {
436: Mat_Pastix *lu = (Mat_Pastix*)(F)->data;
438: lu->iparm[IPARM_FACTORIZATION] = API_FACT_LLT;
439: lu->iparm[IPARM_SYM] = API_SYM_NO;
440: lu->matstruc = DIFFERENT_NONZERO_PATTERN;
441: (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
442: return 0;
443: }
445: PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
446: {
447: PetscBool iascii;
448: PetscViewerFormat format;
450: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
451: if (iascii) {
452: PetscViewerGetFormat(viewer,&format);
453: if (format == PETSC_VIEWER_ASCII_INFO) {
454: Mat_Pastix *lu=(Mat_Pastix*)A->data;
456: PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");
457: PetscViewerASCIIPrintf(viewer," Matrix type : %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
458: PetscViewerASCIIPrintf(viewer," Level of printing (0,1,2): %" PetscInt_FMT " \n",lu->iparm[IPARM_VERBOSE]);
459: PetscViewerASCIIPrintf(viewer," Number of refinements iterations : %" PetscInt_FMT " \n",lu->iparm[IPARM_NBITER]);
460: PetscPrintf(PETSC_COMM_SELF," Error : %g \n",lu->dparm[DPARM_RELATIVE_ERROR]);
461: }
462: }
463: return 0;
464: }
466: /*MC
467: MATSOLVERPASTIX - A solver package providing direct solvers (LU) for distributed
468: and sequential matrices via the external package PaStiX.
470: Use ./configure --download-pastix --download-ptscotch to have PETSc installed with PasTiX
472: Use -pc_type lu -pc_factor_mat_solver_type pastix to use this direct solver
474: Options Database Keys:
475: + -mat_pastix_verbose <0,1,2> - print level
476: - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task.
478: Notes:
479: This only works for matrices with symmetric nonzero structure, if you pass it a matrix with
480: nonsymmetric structure PasTiX and hence PETSc return with an error.
482: Level: beginner
484: .seealso: PCFactorSetMatSolverType(), MatSolverType
486: M*/
488: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
489: {
490: Mat_Pastix *lu =(Mat_Pastix*)A->data;
492: info->block_size = 1.0;
493: info->nz_allocated = lu->iparm[IPARM_NNZEROS];
494: info->nz_used = lu->iparm[IPARM_NNZEROS];
495: info->nz_unneeded = 0.0;
496: info->assemblies = 0.0;
497: info->mallocs = 0.0;
498: info->memory = 0.0;
499: info->fill_ratio_given = 0;
500: info->fill_ratio_needed = 0;
501: info->factor_mallocs = 0;
502: return 0;
503: }
505: static PetscErrorCode MatFactorGetSolverType_pastix(Mat A,MatSolverType *type)
506: {
507: *type = MATSOLVERPASTIX;
508: return 0;
509: }
511: /*
512: The seq and mpi versions of this function are the same
513: */
514: static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
515: {
516: Mat B;
517: Mat_Pastix *pastix;
520: /* Create the factorization matrix */
521: MatCreate(PetscObjectComm((PetscObject)A),&B);
522: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
523: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
524: MatSetUp(B);
526: B->trivialsymbolic = PETSC_TRUE;
527: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
528: B->ops->view = MatView_PaStiX;
529: B->ops->getinfo = MatGetInfo_PaStiX;
531: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);
533: B->factortype = MAT_FACTOR_LU;
535: /* set solvertype */
536: PetscFree(B->solvertype);
537: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
539: PetscNewLog(B,&pastix);
541: pastix->CleanUpPastix = PETSC_FALSE;
542: pastix->scat_rhs = NULL;
543: pastix->scat_sol = NULL;
544: B->ops->getinfo = MatGetInfo_External;
545: B->ops->destroy = MatDestroy_Pastix;
546: B->data = (void*)pastix;
548: *F = B;
549: return 0;
550: }
552: static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
553: {
554: Mat B;
555: Mat_Pastix *pastix;
558: /* Create the factorization matrix */
559: MatCreate(PetscObjectComm((PetscObject)A),&B);
560: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
561: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
562: MatSetUp(B);
564: B->trivialsymbolic = PETSC_TRUE;
565: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
566: B->ops->view = MatView_PaStiX;
567: B->ops->getinfo = MatGetInfo_PaStiX;
568: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);
570: B->factortype = MAT_FACTOR_LU;
572: /* set solvertype */
573: PetscFree(B->solvertype);
574: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
576: PetscNewLog(B,&pastix);
578: pastix->CleanUpPastix = PETSC_FALSE;
579: pastix->scat_rhs = NULL;
580: pastix->scat_sol = NULL;
581: B->ops->getinfo = MatGetInfo_External;
582: B->ops->destroy = MatDestroy_Pastix;
583: B->data = (void*)pastix;
585: *F = B;
586: return 0;
587: }
589: static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
590: {
591: Mat B;
592: Mat_Pastix *pastix;
595: /* Create the factorization matrix */
596: MatCreate(PetscObjectComm((PetscObject)A),&B);
597: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
598: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
599: MatSetUp(B);
601: B->trivialsymbolic = PETSC_TRUE;
602: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
603: B->ops->view = MatView_PaStiX;
604: B->ops->getinfo = MatGetInfo_PaStiX;
605: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);
607: B->factortype = MAT_FACTOR_CHOLESKY;
609: /* set solvertype */
610: PetscFree(B->solvertype);
611: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
613: PetscNewLog(B,&pastix);
615: pastix->CleanUpPastix = PETSC_FALSE;
616: pastix->scat_rhs = NULL;
617: pastix->scat_sol = NULL;
618: B->ops->getinfo = MatGetInfo_External;
619: B->ops->destroy = MatDestroy_Pastix;
620: B->data = (void*)pastix;
621: *F = B;
622: return 0;
623: }
625: static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
626: {
627: Mat B;
628: Mat_Pastix *pastix;
632: /* Create the factorization matrix */
633: MatCreate(PetscObjectComm((PetscObject)A),&B);
634: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
635: PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
636: MatSetUp(B);
638: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
639: B->ops->view = MatView_PaStiX;
640: B->ops->getinfo = MatGetInfo_PaStiX;
641: B->ops->destroy = MatDestroy_Pastix;
642: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);
644: B->factortype = MAT_FACTOR_CHOLESKY;
646: /* set solvertype */
647: PetscFree(B->solvertype);
648: PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);
650: PetscNewLog(B,&pastix);
652: pastix->CleanUpPastix = PETSC_FALSE;
653: pastix->scat_rhs = NULL;
654: pastix->scat_sol = NULL;
655: B->data = (void*)pastix;
657: *F = B;
658: return 0;
659: }
661: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Pastix(void)
662: {
663: MatSolverTypeRegister(MATSOLVERPASTIX,MATMPIAIJ, MAT_FACTOR_LU,MatGetFactor_mpiaij_pastix);
664: MatSolverTypeRegister(MATSOLVERPASTIX,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_pastix);
665: MatSolverTypeRegister(MATSOLVERPASTIX,MATMPISBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_mpisbaij_pastix);
666: MatSolverTypeRegister(MATSOLVERPASTIX,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_pastix);
667: return 0;
668: }