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: }