Actual source code: pastix.c

petsc-master 2017-07-23
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*);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

207:   }
208:   return(0);
209: }

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

221:   if (lu->CleanUpPastix) {
222:     /* Terminate instance, deallocate memories */
223:     VecScatterDestroy(&lu->scat_rhs);
224:     VecDestroy(&lu->b_seq);
225:     VecScatterDestroy(&lu->scat_sol);

227:     lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN;
228:     lu->iparm[IPARM_END_TASK]  =API_TASK_CLEAN;

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

243:     PetscFree(lu->colptr);
244:     PetscFree(lu->row);
245:     PetscFree(lu->val);
246:     PetscFree(lu->perm);
247:     PetscFree(lu->invp);
248:     MPI_Comm_free(&(lu->pastix_comm));
249:   }
250:   PetscFree(A->data);
251:   return(0);
252: }

254: /*
255:   Gather right-hand-side.
256:   Call for Solve step.
257:   Scatter solution.
258:  */
259: PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
260: {
261:   Mat_Pastix     *lu=(Mat_Pastix*)A->data;
262:   PetscScalar    *array;
263:   Vec            x_seq;

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

285:   /* solve phase */
286:   /*-------------*/
287:   lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
288:   lu->iparm[IPARM_END_TASK]   = API_TASK_REFINE;
289:   lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;

291:   PASTIX_CALL(&(lu->pastix_data),
292:               lu->pastix_comm,
293:               lu->n,
294:               lu->colptr,
295:               lu->row,
296:               (PastixScalar*)lu->val,
297:               lu->perm,
298:               lu->invp,
299:               (PastixScalar*)lu->rhs,
300:               lu->rhsnbr,
301:               lu->iparm,
302:               lu->dparm);

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

306:   if (lu->commSize == 1) {
307:     VecRestoreArray(x,&(lu->rhs));
308:   } else {
309:     VecRestoreArray(x_seq,&(lu->rhs));
310:   }

312:   if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
313:     VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
314:     VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
315:   }
316:   return(0);
317: }

319: /*
320:   Numeric factorisation using PaStiX solver.

322:  */
323: PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
324: {
325:   Mat_Pastix     *lu =(Mat_Pastix*)(F)->data;
326:   Mat            *tseq;
327:   PetscErrorCode 0;
328:   PetscInt       icntl;
329:   PetscInt       M=A->rmap->N;
330:   PetscBool      valOnly,flg, isSym;
331:   IS             is_iden;
332:   Vec            b;
333:   IS             isrow;
334:   PetscBool      isSeqAIJ,isSeqSBAIJ,isMPIAIJ;

337:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
338:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
339:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
340:   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
341:     (F)->ops->solve = MatSolve_PaStiX;

343:     /* Initialize a PASTIX instance */
344:     MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm));
345:     MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
346:     MPI_Comm_size(lu->pastix_comm, &lu->commSize);

348:     /* Set pastix options */
349:     lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
350:     lu->iparm[IPARM_START_TASK]       = API_TASK_INIT;
351:     lu->iparm[IPARM_END_TASK]         = API_TASK_INIT;

353:     lu->rhsnbr = 1;

355:     /* Call to set default pastix options */
356:     PASTIX_CALL(&(lu->pastix_data),
357:                 lu->pastix_comm,
358:                 lu->n,
359:                 lu->colptr,
360:                 lu->row,
361:                 (PastixScalar*)lu->val,
362:                 lu->perm,
363:                 lu->invp,
364:                 (PastixScalar*)lu->rhs,
365:                 lu->rhsnbr,
366:                 lu->iparm,
367:                 lu->dparm);

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

371:     icntl = -1;

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

375:     PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg);
376:     if ((flg && icntl >= 0) || PetscLogPrintInfo) {
377:       lu->iparm[IPARM_VERBOSE] =  icntl;
378:     }
379:     icntl=-1;
380:     PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg);
381:     if ((flg && icntl > 0)) {
382:       lu->iparm[IPARM_THREAD_NBR] = icntl;
383:     }
384:     PetscOptionsEnd();
385:     valOnly = PETSC_FALSE;
386:   } else {
387:     if (isSeqAIJ || isMPIAIJ) {
388:       PetscFree(lu->colptr);
389:       PetscFree(lu->row);
390:       PetscFree(lu->val);
391:       valOnly = PETSC_FALSE;
392:     } else valOnly = PETSC_TRUE;
393:   }

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

397:   /* convert mpi A to seq mat A */
398:   ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);
399:   MatCreateSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);
400:   ISDestroy(&isrow);

402:   MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
403:   MatIsSymmetric(*tseq,0.0,&isSym);
404:   MatDestroyMatrices(1,&tseq);

406:   if (!lu->perm) {
407:     PetscMalloc1(lu->n,&(lu->perm));
408:     PetscMalloc1(lu->n,&(lu->invp));
409:   }

411:   if (isSym) {
412:     /* On symmetric matrix, LLT */
413:     lu->iparm[IPARM_SYM]           = API_SYM_YES;
414:     lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
415:   } else {
416:     /* On unsymmetric matrix, LU */
417:     lu->iparm[IPARM_SYM]           = API_SYM_NO;
418:     lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
419:   }

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

436:     PASTIX_CALL(&(lu->pastix_data),
437:                 lu->pastix_comm,
438:                 lu->n,
439:                 lu->colptr,
440:                 lu->row,
441:                 (PastixScalar*)lu->val,
442:                 lu->perm,
443:                 lu->invp,
444:                 (PastixScalar*)lu->rhs,
445:                 lu->rhsnbr,
446:                 lu->iparm,
447:                 lu->dparm);
448:     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]);
449:   } else {
450:     lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
451:     lu->iparm[IPARM_END_TASK]   = API_TASK_NUMFACT;
452:     PASTIX_CALL(&(lu->pastix_data),
453:                 lu->pastix_comm,
454:                 lu->n,
455:                 lu->colptr,
456:                 lu->row,
457:                 (PastixScalar*)lu->val,
458:                 lu->perm,
459:                 lu->invp,
460:                 (PastixScalar*)lu->rhs,
461:                 lu->rhsnbr,
462:                 lu->iparm,
463:                 lu->dparm);
464:     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]);
465:   }

467:   (F)->assembled    = PETSC_TRUE;
468:   lu->matstruc      = SAME_NONZERO_PATTERN;
469:   lu->CleanUpPastix = PETSC_TRUE;
470:   return(0);
471: }

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

479:   lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
480:   lu->iparm[IPARM_SYM]           = API_SYM_YES;
481:   lu->matstruc                   = DIFFERENT_NONZERO_PATTERN;
482:   F->ops->lufactornumeric        = MatFactorNumeric_PaStiX;
483:   return(0);
484: }


487: /* Note the Petsc r permutation is ignored */
488: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
489: {
490:   Mat_Pastix *lu = (Mat_Pastix*)(F)->data;

493:   lu->iparm[IPARM_FACTORIZATION]  = API_FACT_LLT;
494:   lu->iparm[IPARM_SYM]            = API_SYM_NO;
495:   lu->matstruc                    = DIFFERENT_NONZERO_PATTERN;
496:   (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
497:   return(0);
498: }

500: PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
501: {
502:   PetscErrorCode    ierr;
503:   PetscBool         iascii;
504:   PetscViewerFormat format;

507:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
508:   if (iascii) {
509:     PetscViewerGetFormat(viewer,&format);
510:     if (format == PETSC_VIEWER_ASCII_INFO) {
511:       Mat_Pastix *lu=(Mat_Pastix*)A->data;

513:       PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");
514:       PetscViewerASCIIPrintf(viewer,"  Matrix type :                      %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
515:       PetscViewerASCIIPrintf(viewer,"  Level of printing (0,1,2):         %d \n",lu->iparm[IPARM_VERBOSE]);
516:       PetscViewerASCIIPrintf(viewer,"  Number of refinements iterations : %d \n",lu->iparm[IPARM_NBITER]);
517:       PetscPrintf(PETSC_COMM_SELF,"  Error :                        %g \n",lu->dparm[DPARM_RELATIVE_ERROR]);
518:     }
519:   }
520:   return(0);
521: }


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

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

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

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

536:   Level: beginner

538: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage

540: M*/


543: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
544: {
545:   Mat_Pastix *lu =(Mat_Pastix*)A->data;

548:   info->block_size        = 1.0;
549:   info->nz_allocated      = lu->iparm[IPARM_NNZEROS];
550:   info->nz_used           = lu->iparm[IPARM_NNZEROS];
551:   info->nz_unneeded       = 0.0;
552:   info->assemblies        = 0.0;
553:   info->mallocs           = 0.0;
554:   info->memory            = 0.0;
555:   info->fill_ratio_given  = 0;
556:   info->fill_ratio_needed = 0;
557:   info->factor_mallocs    = 0;
558:   return(0);
559: }

561: static PetscErrorCode MatFactorGetSolverPackage_pastix(Mat A,const MatSolverPackage *type)
562: {
564:   *type = MATSOLVERPASTIX;
565:   return(0);
566: }

568: /*
569:     The seq and mpi versions of this function are the same
570: */
571: static PetscErrorCode MatGetFactor_seqaij_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;

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

591:   B->factortype = MAT_FACTOR_LU;
592: 
593:   /* set solvertype */
594:   PetscFree(B->solvertype);
595:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

597:   PetscNewLog(B,&pastix);

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

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

610: static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
611: {
612:   Mat            B;
614:   Mat_Pastix     *pastix;

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

624:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
625:   B->ops->view             = MatView_PaStiX;
626:   B->ops->getinfo          = MatGetInfo_PaStiX;
627:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

629:   B->factortype = MAT_FACTOR_LU;

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

635:   PetscNewLog(B,&pastix);

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

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

648: static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
649: {
650:   Mat            B;
652:   Mat_Pastix     *pastix;

655:   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:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

667:   B->factortype = MAT_FACTOR_CHOLESKY;

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

673:   PetscNewLog(B,&pastix);

675:   pastix->CleanUpPastix = PETSC_FALSE;
676:   pastix->scat_rhs      = NULL;
677:   pastix->scat_sol      = NULL;
678:   B->ops->getinfo       = MatGetInfo_External;
679:   B->ops->destroy       = MatDestroy_Pastix;
680:   B->data               = (void*)pastix;

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

686: static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
687: {
688:   Mat            B;
690:   Mat_Pastix     *pastix;

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

695:   /* Create the factorization matrix */
696:   MatCreate(PetscObjectComm((PetscObject)A),&B);
697:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
698:   PetscStrallocpy("pastix",&((PetscObject)B)->type_name);
699:   MatSetUp(B);

701:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
702:   B->ops->view                   = MatView_PaStiX;
703:   B->ops->getinfo                = MatGetInfo_PaStiX;
704:   B->ops->destroy                = MatDestroy_Pastix;
705:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_pastix);

707:   B->factortype = MAT_FACTOR_CHOLESKY;

709:   /* set solvertype */
710:   PetscFree(B->solvertype);
711:   PetscStrallocpy(MATSOLVERPASTIX,&B->solvertype);

713:   PetscNewLog(B,&pastix);

715:   pastix->CleanUpPastix = PETSC_FALSE;
716:   pastix->scat_rhs      = NULL;
717:   pastix->scat_sol      = NULL;
718:   B->data               = (void*)pastix;

720:   *F = B;
721:   return(0);
722: }

724: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Pastix(void)
725: {

729:   MatSolverPackageRegister(MATSOLVERPASTIX,MATMPIAIJ,        MAT_FACTOR_LU,MatGetFactor_mpiaij_pastix);
730:   MatSolverPackageRegister(MATSOLVERPASTIX,MATSEQAIJ,        MAT_FACTOR_LU,MatGetFactor_seqaij_pastix);
731:   MatSolverPackageRegister(MATSOLVERPASTIX,MATMPISBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_mpisbaij_pastix);
732:   MatSolverPackageRegister(MATSOLVERPASTIX,MATSEQSBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_pastix);
733:   return(0);
734: }