Actual source code: pastix.c

petsc-master 2016-12-08
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: #define PASTIX_CHECKMATRIX c_pastix_checkMatrix
 21: #else
 22: #define PASTIX_CALL z_pastix
 23: #define PASTIX_CHECKMATRIX z_pastix_checkMatrix
 24: #endif

 26: #else /* PETSC_USE_COMPLEX */

 28: #if defined(PETSC_USE_REAL_SINGLE)
 29: #define PASTIX_CALL s_pastix
 30: #define PASTIX_CHECKMATRIX s_pastix_checkMatrix
 31: #else
 32: #define PASTIX_CALL d_pastix
 33: #define PASTIX_CHECKMATRIX d_pastix_checkMatrix
 34: #endif

 36: #endif /* PETSC_USE_COMPLEX */

 38: typedef PetscScalar PastixScalar;

 40: typedef struct Mat_Pastix_ {
 41:   pastix_data_t *pastix_data;    /* Pastix data storage structure                        */
 42:   MatStructure  matstruc;
 43:   PetscInt      n;               /* Number of columns in the matrix                      */
 44:   PetscInt      *colptr;         /* Index of first element of each column in row and val */
 45:   PetscInt      *row;            /* Row of each element of the matrix                    */
 46:   PetscScalar   *val;            /* Value of each element of the matrix                  */
 47:   PetscInt      *perm;           /* Permutation tabular                                  */
 48:   PetscInt      *invp;           /* Reverse permutation tabular                          */
 49:   PetscScalar   *rhs;            /* Rhight-hand-side member                              */
 50:   PetscInt      rhsnbr;          /* Rhight-hand-side number (must be 1)                  */
 51:   PetscInt      iparm[64];       /* Integer parameters                                   */
 52:   double        dparm[64];       /* Floating point parameters                            */
 53:   MPI_Comm      pastix_comm;     /* PaStiX MPI communicator                              */
 54:   PetscMPIInt   commRank;        /* MPI rank                                             */
 55:   PetscMPIInt   commSize;        /* MPI communicator size                                */
 56:   PetscBool     CleanUpPastix;   /* Boolean indicating if we call PaStiX clean step      */
 57:   VecScatter    scat_rhs;
 58:   VecScatter    scat_sol;
 59:   Vec           b_seq;
 60: } Mat_Pastix;

 62: extern PetscErrorCode MatDuplicate_Pastix(Mat,MatDuplicateOption,Mat*);

 66: /*
 67:    convert Petsc seqaij matrix to CSC: colptr[n], row[nz], val[nz]

 69:   input:
 70:     A       - matrix in seqaij or mpisbaij (bs=1) format
 71:     valOnly - FALSE: spaces are allocated and values are set for the CSC
 72:               TRUE:  Only fill values
 73:   output:
 74:     n       - Size of the matrix
 75:     colptr  - Index of first element of each column in row and val
 76:     row     - Row of each element of the matrix
 77:     values  - Value of each element of the matrix
 78:  */
 79: PetscErrorCode MatConvertToCSC(Mat A,PetscBool valOnly,PetscInt *n,PetscInt **colptr,PetscInt **row,PetscScalar **values)
 80: {
 81:   Mat_SeqAIJ     *aa      = (Mat_SeqAIJ*)A->data;
 82:   PetscInt       *rowptr  = aa->i;
 83:   PetscInt       *col     = aa->j;
 84:   PetscScalar    *rvalues = aa->a;
 85:   PetscInt       m        = A->rmap->N;
 86:   PetscInt       nnz;
 87:   PetscInt       i,j, k;
 88:   PetscInt       base = 1;
 89:   PetscInt       idx;
 91:   PetscInt       colidx;
 92:   PetscInt       *colcount;
 93:   PetscBool      isSBAIJ;
 94:   PetscBool      isSeqSBAIJ;
 95:   PetscBool      isMpiSBAIJ;
 96:   PetscBool      isSym;
 97:   PetscBool      flg;
 98:   PetscInt       icntl;
 99:   PetscInt       verb;
100:   PetscInt       check;

103:   MatIsSymmetric(A,0.0,&isSym);
104:   PetscObjectTypeCompare((PetscObject)A,MATSBAIJ,&isSBAIJ);
105:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
106:   PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMpiSBAIJ);

108:   *n = A->cmap->N;

110:   /* PaStiX only needs triangular matrix if matrix is symmetric
111:    */
112:   if (isSym && !(isSBAIJ || isSeqSBAIJ || isMpiSBAIJ)) nnz = (aa->nz - *n)/2 + *n;
113:   else nnz = aa->nz;

115:   if (!valOnly) {
116:     PetscMalloc1((*n)+1,colptr);
117:     PetscMalloc1(nnz,row);
118:     PetscMalloc1(nnz,values);

120:     if (isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) {
121:       PetscMemcpy (*colptr, rowptr, ((*n)+1)*sizeof(PetscInt));
122:       for (i = 0; i < *n+1; i++) (*colptr)[i] += base;
123:       PetscMemcpy (*row, col, (nnz)*sizeof(PetscInt));
124:       for (i = 0; i < nnz; i++) (*row)[i] += base;
125:       PetscMemcpy (*values, rvalues, (nnz)*sizeof(PetscScalar));
126:     } else {
127:       PetscMalloc1(*n,&colcount);

129:       for (i = 0; i < m; i++) colcount[i] = 0;
130:       /* Fill-in colptr */
131:       for (i = 0; i < m; i++) {
132:         for (j = rowptr[i]; j < rowptr[i+1]; j++) {
133:           if (!isSym || col[j] <= i)  colcount[col[j]]++;
134:         }
135:       }

137:       (*colptr)[0] = base;
138:       for (j = 0; j < *n; j++) {
139:         (*colptr)[j+1] = (*colptr)[j] + colcount[j];
140:         /* in next loop we fill starting from (*colptr)[colidx] - base */
141:         colcount[j] = -base;
142:       }

144:       /* Fill-in rows and values */
145:       for (i = 0; i < m; i++) {
146:         for (j = rowptr[i]; j < rowptr[i+1]; j++) {
147:           if (!isSym || col[j] <= i) {
148:             colidx         = col[j];
149:             idx            = (*colptr)[colidx] + colcount[colidx];
150:             (*row)[idx]    = i + base;
151:             (*values)[idx] = rvalues[j];
152:             colcount[colidx]++;
153:           }
154:         }
155:       }
156:       PetscFree(colcount);
157:     }
158:   } else {
159:     /* Fill-in only values */
160:     for (i = 0; i < m; i++) {
161:       for (j = rowptr[i]; j < rowptr[i+1]; j++) {
162:         colidx = col[j];
163:         if ((isSBAIJ || isSeqSBAIJ || isMpiSBAIJ) ||!isSym || col[j] <= i) {
164:           /* look for the value to fill */
165:           for (k = (*colptr)[colidx] - base; k < (*colptr)[colidx + 1] - base; k++) {
166:             if (((*row)[k]-base) == i) {
167:               (*values)[k] = rvalues[j];
168:               break;
169:             }
170:           }
171:           /* data structure of sparse matrix has changed */
172:           if (k == (*colptr)[colidx + 1] - base) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"overflow on k %D",k);
173:         }
174:       }
175:     }
176:   }

178:   icntl =-1;
179:   check = 0;
180:   PetscOptionsGetInt(NULL,((PetscObject) A)->prefix, "-mat_pastix_check", &icntl, &flg);
181:   if ((flg && icntl >= 0) || PetscLogPrintInfo) check =  icntl;

183:   if (check == 1) {
184:     PetscScalar *tmpvalues;
185:     PetscInt    *tmprows,*tmpcolptr;

187:     PetscMalloc3(nnz,&tmpvalues,nnz,&tmprows,*n+1,&tmpcolptr);

189:     PetscMemcpy(tmpcolptr,*colptr,(*n+1)*sizeof(PetscInt));
190:     PetscMemcpy(tmprows,*row,nnz*sizeof(PetscInt));
191:     PetscMemcpy(tmpvalues,*values,nnz*sizeof(PetscScalar));
192:     PetscFree(*row);
193:     PetscFree(*values);

195:     icntl=-1;
196:     verb = API_VERBOSE_NOT;
197:     /* "iparm[IPARM_VERBOSE] : level of printing (0 to 2)" */
198:     PetscOptionsGetInt(NULL,((PetscObject) A)->prefix, "-mat_pastix_verbose", &icntl, &flg);
199:     if ((flg && icntl >= 0) || PetscLogPrintInfo) verb =  icntl;
200:     PASTIX_CHECKMATRIX(MPI_COMM_WORLD,verb,((isSym != 0) ? API_SYM_YES : API_SYM_NO),API_YES,*n,&tmpcolptr,&tmprows,(PastixScalar**)&tmpvalues,NULL,1);

202:     PetscMemcpy(*colptr,tmpcolptr,(*n+1)*sizeof(PetscInt));
203:     PetscMalloc1(((*colptr)[*n]-1),row);
204:     PetscMemcpy(*row,tmprows,((*colptr)[*n]-1)*sizeof(PetscInt));
205:     PetscMalloc1(((*colptr)[*n]-1),values);
206:     PetscMemcpy(*values,tmpvalues,((*colptr)[*n]-1)*sizeof(PetscScalar));
207:     PetscFree3(tmpvalues,tmprows,tmpcolptr);

209:   }
210:   return(0);
211: }

215: /*
216:   Call clean step of PaStiX if lu->CleanUpPastix == true.
217:   Free the CSC matrix.
218:  */
219: PetscErrorCode MatDestroy_Pastix(Mat A)
220: {
221:   Mat_Pastix     *lu=(Mat_Pastix*)A->data;

225:   if (lu->CleanUpPastix) {
226:     /* Terminate instance, deallocate memories */
227:     VecScatterDestroy(&lu->scat_rhs);
228:     VecDestroy(&lu->b_seq);
229:     VecScatterDestroy(&lu->scat_sol);

231:     lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN;
232:     lu->iparm[IPARM_END_TASK]  =API_TASK_CLEAN;

234:     PASTIX_CALL(&(lu->pastix_data),
235:                 lu->pastix_comm,
236:                 lu->n,
237:                 lu->colptr,
238:                 lu->row,
239:                 (PastixScalar*)lu->val,
240:                 lu->perm,
241:                 lu->invp,
242:                 (PastixScalar*)lu->rhs,
243:                 lu->rhsnbr,
244:                 lu->iparm,
245:                 lu->dparm);

247:     PetscFree(lu->colptr);
248:     PetscFree(lu->row);
249:     PetscFree(lu->val);
250:     PetscFree(lu->perm);
251:     PetscFree(lu->invp);
252:     MPI_Comm_free(&(lu->pastix_comm));
253:   }
254:   PetscFree(A->data);
255:   return(0);
256: }

260: /*
261:   Gather right-hand-side.
262:   Call for Solve step.
263:   Scatter solution.
264:  */
265: PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
266: {
267:   Mat_Pastix     *lu=(Mat_Pastix*)A->data;
268:   PetscScalar    *array;
269:   Vec            x_seq;

273:   lu->rhsnbr = 1;
274:   x_seq      = lu->b_seq;
275:   if (lu->commSize > 1) {
276:     /* PaStiX only supports centralized rhs. Scatter b into a seqential rhs vector */
277:     VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
278:     VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
279:     VecGetArray(x_seq,&array);
280:   } else {  /* size == 1 */
281:     VecCopy(b,x);
282:     VecGetArray(x,&array);
283:   }
284:   lu->rhs = array;
285:   if (lu->commSize == 1) {
286:     VecRestoreArray(x,&array);
287:   } else {
288:     VecRestoreArray(x_seq,&array);
289:   }

291:   /* solve phase */
292:   /*-------------*/
293:   lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
294:   lu->iparm[IPARM_END_TASK]   = API_TASK_REFINE;
295:   lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;

297:   PASTIX_CALL(&(lu->pastix_data),
298:               lu->pastix_comm,
299:               lu->n,
300:               lu->colptr,
301:               lu->row,
302:               (PastixScalar*)lu->val,
303:               lu->perm,
304:               lu->invp,
305:               (PastixScalar*)lu->rhs,
306:               lu->rhsnbr,
307:               lu->iparm,
308:               lu->dparm);

310:   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]);

312:   if (lu->commSize == 1) {
313:     VecRestoreArray(x,&(lu->rhs));
314:   } else {
315:     VecRestoreArray(x_seq,&(lu->rhs));
316:   }

318:   if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
319:     VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
320:     VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
321:   }
322:   return(0);
323: }

325: /*
326:   Numeric factorisation using PaStiX solver.

328:  */
331: PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
332: {
333:   Mat_Pastix     *lu =(Mat_Pastix*)(F)->data;
334:   Mat            *tseq;
335:   PetscErrorCode 0;
336:   PetscInt       icntl;
337:   PetscInt       M=A->rmap->N;
338:   PetscBool      valOnly,flg, isSym;
339:   IS             is_iden;
340:   Vec            b;
341:   IS             isrow;
342:   PetscBool      isSeqAIJ,isSeqSBAIJ,isMPIAIJ;

345:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
346:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
347:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
348:   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
349:     (F)->ops->solve = MatSolve_PaStiX;

351:     /* Initialize a PASTIX instance */
352:     MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm));
353:     MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
354:     MPI_Comm_size(lu->pastix_comm, &lu->commSize);

356:     /* Set pastix options */
357:     lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
358:     lu->iparm[IPARM_START_TASK]       = API_TASK_INIT;
359:     lu->iparm[IPARM_END_TASK]         = API_TASK_INIT;

361:     lu->rhsnbr = 1;

363:     /* Call to set default pastix options */
364:     PASTIX_CALL(&(lu->pastix_data),
365:                 lu->pastix_comm,
366:                 lu->n,
367:                 lu->colptr,
368:                 lu->row,
369:                 (PastixScalar*)lu->val,
370:                 lu->perm,
371:                 lu->invp,
372:                 (PastixScalar*)lu->rhs,
373:                 lu->rhsnbr,
374:                 lu->iparm,
375:                 lu->dparm);

377:     PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"PaStiX Options","Mat");

379:     icntl = -1;

381:     lu->iparm[IPARM_VERBOSE] = API_VERBOSE_NOT;

383:     PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg);
384:     if ((flg && icntl >= 0) || PetscLogPrintInfo) {
385:       lu->iparm[IPARM_VERBOSE] =  icntl;
386:     }
387:     icntl=-1;
388:     PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg);
389:     if ((flg && icntl > 0)) {
390:       lu->iparm[IPARM_THREAD_NBR] = icntl;
391:     }
392:     PetscOptionsEnd();
393:     valOnly = PETSC_FALSE;
394:   } else {
395:     if (isSeqAIJ || isMPIAIJ) {
396:       PetscFree(lu->colptr);
397:       PetscFree(lu->row);
398:       PetscFree(lu->val);
399:       valOnly = PETSC_FALSE;
400:     } else valOnly = PETSC_TRUE;
401:   }

403:   lu->iparm[IPARM_MATRIX_VERIFICATION] = API_YES;

405:   /* convert mpi A to seq mat A */
406:   ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);
407:   MatGetSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);
408:   ISDestroy(&isrow);

410:   MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
411:   MatIsSymmetric(*tseq,0.0,&isSym);
412:   MatDestroyMatrices(1,&tseq);

414:   if (!lu->perm) {
415:     PetscMalloc1(lu->n,&(lu->perm));
416:     PetscMalloc1(lu->n,&(lu->invp));
417:   }

419:   if (isSym) {
420:     /* On symmetric matrix, LLT */
421:     lu->iparm[IPARM_SYM]           = API_SYM_YES;
422:     lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
423:   } else {
424:     /* On unsymmetric matrix, LU */
425:     lu->iparm[IPARM_SYM]           = API_SYM_NO;
426:     lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
427:   }

429:   /*----------------*/
430:   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
431:     if (!(isSeqAIJ || isSeqSBAIJ) && !lu->b_seq) {
432:       /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
433:       VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);
434:       ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
435:       MatCreateVecs(A,NULL,&b);
436:       VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);
437:       VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol);
438:       ISDestroy(&is_iden);
439:       VecDestroy(&b);
440:     }
441:     lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
442:     lu->iparm[IPARM_END_TASK]   = API_TASK_NUMFACT;

444:     PASTIX_CALL(&(lu->pastix_data),
445:                 lu->pastix_comm,
446:                 lu->n,
447:                 lu->colptr,
448:                 lu->row,
449:                 (PastixScalar*)lu->val,
450:                 lu->perm,
451:                 lu->invp,
452:                 (PastixScalar*)lu->rhs,
453:                 lu->rhsnbr,
454:                 lu->iparm,
455:                 lu->dparm);
456:     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]);
457:   } else {
458:     lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
459:     lu->iparm[IPARM_END_TASK]   = API_TASK_NUMFACT;
460:     PASTIX_CALL(&(lu->pastix_data),
461:                 lu->pastix_comm,
462:                 lu->n,
463:                 lu->colptr,
464:                 lu->row,
465:                 (PastixScalar*)lu->val,
466:                 lu->perm,
467:                 lu->invp,
468:                 (PastixScalar*)lu->rhs,
469:                 lu->rhsnbr,
470:                 lu->iparm,
471:                 lu->dparm);
472:     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]);
473:   }

475:   (F)->assembled    = PETSC_TRUE;
476:   lu->matstruc      = SAME_NONZERO_PATTERN;
477:   lu->CleanUpPastix = PETSC_TRUE;
478:   return(0);
479: }

481: /* Note the Petsc r and c permutations are ignored */
484: PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
485: {
486:   Mat_Pastix *lu = (Mat_Pastix*)F->data;

489:   lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
490:   lu->iparm[IPARM_SYM]           = API_SYM_YES;
491:   lu->matstruc                   = DIFFERENT_NONZERO_PATTERN;
492:   F->ops->lufactornumeric        = MatFactorNumeric_PaStiX;
493:   return(0);
494: }


497: /* Note the Petsc r permutation is ignored */
500: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
501: {
502:   Mat_Pastix *lu = (Mat_Pastix*)(F)->data;

505:   lu->iparm[IPARM_FACTORIZATION]  = API_FACT_LLT;
506:   lu->iparm[IPARM_SYM]            = API_SYM_NO;
507:   lu->matstruc                    = DIFFERENT_NONZERO_PATTERN;
508:   (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
509:   return(0);
510: }

514: PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
515: {
516:   PetscErrorCode    ierr;
517:   PetscBool         iascii;
518:   PetscViewerFormat format;

521:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
522:   if (iascii) {
523:     PetscViewerGetFormat(viewer,&format);
524:     if (format == PETSC_VIEWER_ASCII_INFO) {
525:       Mat_Pastix *lu=(Mat_Pastix*)A->data;

527:       PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");
528:       PetscViewerASCIIPrintf(viewer,"  Matrix type :                      %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
529:       PetscViewerASCIIPrintf(viewer,"  Level of printing (0,1,2):         %d \n",lu->iparm[IPARM_VERBOSE]);
530:       PetscViewerASCIIPrintf(viewer,"  Number of refinements iterations : %d \n",lu->iparm[IPARM_NBITER]);
531:       PetscPrintf(PETSC_COMM_SELF,"  Error :                        %g \n",lu->dparm[DPARM_RELATIVE_ERROR]);
532:     }
533:   }
534:   return(0);
535: }


538: /*MC
539:      MATSOLVERPASTIX  - A solver package providing direct solvers (LU) for distributed
540:   and sequential matrices via the external package PaStiX.

542:   Use ./configure --download-pastix --download-ptscotch  to have PETSc installed with PasTiX

544:   Use -pc_type lu -pc_factor_mat_solver_package pastix to us this direct solver

546:   Options Database Keys:
547: + -mat_pastix_verbose   <0,1,2>   - print level
548: - -mat_pastix_threadnbr <integer> - Set the thread number by MPI task.

550:   Level: beginner

552: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage

554: M*/


559: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
560: {
561:   Mat_Pastix *lu =(Mat_Pastix*)A->data;

564:   info->block_size        = 1.0;
565:   info->nz_allocated      = lu->iparm[IPARM_NNZEROS];
566:   info->nz_used           = lu->iparm[IPARM_NNZEROS];
567:   info->nz_unneeded       = 0.0;
568:   info->assemblies        = 0.0;
569:   info->mallocs           = 0.0;
570:   info->memory            = 0.0;
571:   info->fill_ratio_given  = 0;
572:   info->fill_ratio_needed = 0;
573:   info->factor_mallocs    = 0;
574:   return(0);
575: }

579: static PetscErrorCode MatFactorGetSolverPackage_pastix(Mat A,const MatSolverPackage *type)
580: {
582:   *type = MATSOLVERPASTIX;
583:   return(0);
584: }

586: /*
587:     The seq and mpi versions of this function are the same
588: */
591: static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
592: {
593:   Mat            B;
595:   Mat_Pastix     *pastix;

598:   if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
599:   /* Create the factorization matrix */
600:   MatCreate(PetscObjectComm((PetscObject)A),&B);
601:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
602:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
603:   MatSetUp(B);

605:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
606:   B->ops->view             = MatView_PaStiX;
607:   B->ops->getinfo          = MatGetInfo_PaStiX;

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

611:   B->factortype = MAT_FACTOR_LU;
612: 
613:   /* set solvertype */
614:   PetscFree(B->solvertype);
615:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

617:   PetscNewLog(B,&pastix);

619:   pastix->CleanUpPastix = PETSC_FALSE;
620:   pastix->scat_rhs      = NULL;
621:   pastix->scat_sol      = NULL;
622:   B->ops->getinfo       = MatGetInfo_External;
623:   B->ops->destroy       = MatDestroy_Pastix;
624:   B->data               = (void*)pastix;

626:   *F = B;
627:   return(0);
628: }

632: static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
633: {
634:   Mat            B;
636:   Mat_Pastix     *pastix;

639:   if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
640:   /* Create the factorization matrix */
641:   MatCreate(PetscObjectComm((PetscObject)A),&B);
642:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
643:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
644:   MatSetUp(B);

646:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
647:   B->ops->view             = MatView_PaStiX;
648:   B->ops->getinfo          = MatGetInfo_PaStiX;
649:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

651:   B->factortype = MAT_FACTOR_LU;

653:   /* set solvertype */
654:   PetscFree(B->solvertype);
655:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

657:   PetscNewLog(B,&pastix);

659:   pastix->CleanUpPastix = PETSC_FALSE;
660:   pastix->scat_rhs      = NULL;
661:   pastix->scat_sol      = NULL;
662:   B->ops->getinfo       = MatGetInfo_External;
663:   B->ops->destroy       = MatDestroy_Pastix;
664:   B->data               = (void*)pastix;

666:   *F = B;
667:   return(0);
668: }

672: static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
673: {
674:   Mat            B;
676:   Mat_Pastix     *pastix;

679:   if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
680:   /* Create the factorization matrix */
681:   MatCreate(PetscObjectComm((PetscObject)A),&B);
682:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
683:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
684:   MatSetUp(B);

686:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
687:   B->ops->view                   = MatView_PaStiX;
688:   B->ops->getinfo                = MatGetInfo_PaStiX;
689:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

691:   B->factortype = MAT_FACTOR_CHOLESKY;

693:   /* set solvertype */
694:   PetscFree(B->solvertype);
695:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

697:   PetscNewLog(B,&pastix);

699:   pastix->CleanUpPastix = PETSC_FALSE;
700:   pastix->scat_rhs      = NULL;
701:   pastix->scat_sol      = NULL;
702:   B->ops->getinfo       = MatGetInfo_External;
703:   B->ops->destroy       = MatDestroy_Pastix;
704:   B->data               = (void*)pastix;

706:   *F = B;
707:   return(0);
708: }

712: static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
713: {
714:   Mat            B;
716:   Mat_Pastix     *pastix;

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

721:   /* Create the factorization matrix */
722:   MatCreate(PetscObjectComm((PetscObject)A),&B);
723:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
724:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
725:   MatSetUp(B);

727:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
728:   B->ops->view                   = MatView_PaStiX;
729:   B->ops->getinfo                = MatGetInfo_PaStiX;
730:   B->ops->destroy                = MatDestroy_Pastix;
731:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

733:   B->factortype = MAT_FACTOR_CHOLESKY;

735:   /* set solvertype */
736:   PetscFree(B->solvertype);
737:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

739:   PetscNewLog(B,&pastix);

741:   pastix->CleanUpPastix = PETSC_FALSE;
742:   pastix->scat_rhs      = NULL;
743:   pastix->scat_sol      = NULL;
744:   B->data               = (void*)pastix;

746:   *F = B;
747:   return(0);
748: }

752: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Pastix(void)
753: {

757:   MatSolverPackageRegister(MATSOLVERPASTIX,MATMPIAIJ,        MAT_FACTOR_LU,MatGetFactor_mpiaij_pastix);
758:   MatSolverPackageRegister(MATSOLVERPASTIX,MATSEQAIJ,        MAT_FACTOR_LU,MatGetFactor_seqaij_pastix);
759:   MatSolverPackageRegister(MATSOLVERPASTIX,MATMPISBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_mpisbaij_pastix);
760:   MatSolverPackageRegister(MATSOLVERPASTIX,MATSEQSBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_pastix);
761:   return(0);
762: }