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;
 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:       PetscArraycpy (*colptr, rowptr, (*n)+1);
112:       for (i = 0; i < *n+1; i++) (*colptr)[i] += base;
113:       PetscArraycpy (*row, col, nnz);
114:       for (i = 0; i < nnz; i++) (*row)[i] += base;
115:       PetscArraycpy (*values, rvalues, nnz);
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:   B->useordering = PETSC_TRUE;

554:   /* set solvertype */
555:   PetscFree(B->solvertype);
556:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

558:   PetscNewLog(B,&pastix);

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

567:   *F = B;
568:   return(0);
569: }

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

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

585:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
586:   B->ops->view             = MatView_PaStiX;
587:   B->ops->getinfo          = MatGetInfo_PaStiX;
588:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);

590:   B->factortype = MAT_FACTOR_LU;

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

596:   PetscNewLog(B,&pastix);

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

605:   *F = B;
606:   return(0);
607: }

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

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

623:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
624:   B->ops->view                   = MatView_PaStiX;
625:   B->ops->getinfo                = MatGetInfo_PaStiX;
626:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);

628:   B->factortype = MAT_FACTOR_CHOLESKY;

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

634:   PetscNewLog(B,&pastix);

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

643:   *F = B;
644:   return(0);
645: }

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

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

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

662:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
663:   B->ops->view                   = MatView_PaStiX;
664:   B->ops->getinfo                = MatGetInfo_PaStiX;
665:   B->ops->destroy                = MatDestroy_Pastix;
666:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_pastix);

668:   B->factortype = MAT_FACTOR_CHOLESKY;

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

674:   PetscNewLog(B,&pastix);

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

681:   *F = B;
682:   return(0);
683: }

685: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Pastix(void)
686: {

690:   MatSolverTypeRegister(MATSOLVERPASTIX,MATMPIAIJ,        MAT_FACTOR_LU,MatGetFactor_mpiaij_pastix);
691:   MatSolverTypeRegister(MATSOLVERPASTIX,MATSEQAIJ,        MAT_FACTOR_LU,MatGetFactor_seqaij_pastix);
692:   MatSolverTypeRegister(MATSOLVERPASTIX,MATMPISBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_mpisbaij_pastix);
693:   MatSolverTypeRegister(MATSOLVERPASTIX,MATSEQSBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_pastix);
694:   return(0);
695: }