Actual source code: pastix.c

petsc-master 2018-01-18
Report Typos and Errors
  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_package pastix to us 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: This only works for matrices with symmetric nonzero structure, if you pass it a matrix with
493:    nonsymmetric structure PasTiX and hence PETSc return with an error.

495:   Level: beginner

497: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage

499: M*/


502: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
503: {
504:   Mat_Pastix *lu =(Mat_Pastix*)A->data;

507:   info->block_size        = 1.0;
508:   info->nz_allocated      = lu->iparm[IPARM_NNZEROS];
509:   info->nz_used           = lu->iparm[IPARM_NNZEROS];
510:   info->nz_unneeded       = 0.0;
511:   info->assemblies        = 0.0;
512:   info->mallocs           = 0.0;
513:   info->memory            = 0.0;
514:   info->fill_ratio_given  = 0;
515:   info->fill_ratio_needed = 0;
516:   info->factor_mallocs    = 0;
517:   return(0);
518: }

520: static PetscErrorCode MatFactorGetSolverPackage_pastix(Mat A,const MatSolverPackage *type)
521: {
523:   *type = MATSOLVERPASTIX;
524:   return(0);
525: }

527: /*
528:     The seq and mpi versions of this function are the same
529: */
530: static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
531: {
532:   Mat            B;
534:   Mat_Pastix     *pastix;

537:   if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
538:   /* Create the factorization matrix */
539:   MatCreate(PetscObjectComm((PetscObject)A),&B);
540:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
541:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
542:   MatSetUp(B);

544:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
545:   B->ops->view             = MatView_PaStiX;
546:   B->ops->getinfo          = MatGetInfo_PaStiX;

548:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

550:   B->factortype = MAT_FACTOR_LU;
551: 
552:   /* set solvertype */
553:   PetscFree(B->solvertype);
554:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

556:   PetscNewLog(B,&pastix);

558:   pastix->CleanUpPastix = PETSC_FALSE;
559:   pastix->scat_rhs      = NULL;
560:   pastix->scat_sol      = NULL;
561:   B->ops->getinfo       = MatGetInfo_External;
562:   B->ops->destroy       = MatDestroy_Pastix;
563:   B->data               = (void*)pastix;

565:   *F = B;
566:   return(0);
567: }

569: static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
570: {
571:   Mat            B;
573:   Mat_Pastix     *pastix;

576:   if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
577:   /* Create the factorization matrix */
578:   MatCreate(PetscObjectComm((PetscObject)A),&B);
579:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
580:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
581:   MatSetUp(B);

583:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
584:   B->ops->view             = MatView_PaStiX;
585:   B->ops->getinfo          = MatGetInfo_PaStiX;
586:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

588:   B->factortype = MAT_FACTOR_LU;

590:   /* set solvertype */
591:   PetscFree(B->solvertype);
592:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

594:   PetscNewLog(B,&pastix);

596:   pastix->CleanUpPastix = PETSC_FALSE;
597:   pastix->scat_rhs      = NULL;
598:   pastix->scat_sol      = NULL;
599:   B->ops->getinfo       = MatGetInfo_External;
600:   B->ops->destroy       = MatDestroy_Pastix;
601:   B->data               = (void*)pastix;

603:   *F = B;
604:   return(0);
605: }

607: static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
608: {
609:   Mat            B;
611:   Mat_Pastix     *pastix;

614:   if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
615:   /* Create the factorization matrix */
616:   MatCreate(PetscObjectComm((PetscObject)A),&B);
617:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
618:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
619:   MatSetUp(B);

621:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
622:   B->ops->view                   = MatView_PaStiX;
623:   B->ops->getinfo                = MatGetInfo_PaStiX;
624:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

626:   B->factortype = MAT_FACTOR_CHOLESKY;

628:   /* set solvertype */
629:   PetscFree(B->solvertype);
630:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

632:   PetscNewLog(B,&pastix);

634:   pastix->CleanUpPastix = PETSC_FALSE;
635:   pastix->scat_rhs      = NULL;
636:   pastix->scat_sol      = NULL;
637:   B->ops->getinfo       = MatGetInfo_External;
638:   B->ops->destroy       = MatDestroy_Pastix;
639:   B->data               = (void*)pastix;

641:   *F = B;
642:   return(0);
643: }

645: static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
646: {
647:   Mat            B;
649:   Mat_Pastix     *pastix;

652:   if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");

654:   /* Create the factorization matrix */
655:   MatCreate(PetscObjectComm((PetscObject)A),&B);
656:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
657:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
658:   MatSetUp(B);

660:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
661:   B->ops->view                   = MatView_PaStiX;
662:   B->ops->getinfo                = MatGetInfo_PaStiX;
663:   B->ops->destroy                = MatDestroy_Pastix;
664:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

666:   B->factortype = MAT_FACTOR_CHOLESKY;

668:   /* set solvertype */
669:   PetscFree(B->solvertype);
670:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

672:   PetscNewLog(B,&pastix);

674:   pastix->CleanUpPastix = PETSC_FALSE;
675:   pastix->scat_rhs      = NULL;
676:   pastix->scat_sol      = NULL;
677:   B->data               = (void*)pastix;

679:   *F = B;
680:   return(0);
681: }

683: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Pastix(void)
684: {

688:   MatSolverPackageRegister(MATSOLVERPASTIX,MATMPIAIJ,        MAT_FACTOR_LU,MatGetFactor_mpiaij_pastix);
689:   MatSolverPackageRegister(MATSOLVERPASTIX,MATSEQAIJ,        MAT_FACTOR_LU,MatGetFactor_seqaij_pastix);
690:   MatSolverPackageRegister(MATSOLVERPASTIX,MATMPISBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_mpisbaij_pastix);
691:   MatSolverPackageRegister(MATSOLVERPASTIX,MATSEQSBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_pastix);
692:   return(0);
693: }