Actual source code: pastix.c

petsc-3.5.2 2014-09-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:   PetscBool     isAIJ;
 61:   PetscErrorCode (*Destroy)(Mat);
 62: } Mat_Pastix;

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

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

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

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

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

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

117:   if (!valOnly) {
118:     PetscMalloc(((*n)+1) *sizeof(PetscInt)   ,colptr);
119:     PetscMalloc(nnz      *sizeof(PetscInt)   ,row);
120:     PetscMalloc1(nnz      ,values);

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

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

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

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

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

185:   if (check == 1) {
186:     PetscScalar *tmpvalues;
187:     PetscInt    *tmprows,*tmpcolptr;
188:     tmpvalues = (PetscScalar*)malloc(nnz*sizeof(PetscScalar));    if (!tmpvalues) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");
189:     tmprows   = (PetscInt*)   malloc(nnz*sizeof(PetscInt));       if (!tmprows)   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");
190:     tmpcolptr = (PetscInt*)   malloc((*n+1)*sizeof(PetscInt));    if (!tmpcolptr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Unable to allocate memory");

192:     PetscMemcpy(tmpcolptr,*colptr,(*n+1)*sizeof(PetscInt));
193:     PetscMemcpy(tmprows,*row,nnz*sizeof(PetscInt));
194:     PetscMemcpy(tmpvalues,*values,nnz*sizeof(PetscScalar));
195:     PetscFree(*row);
196:     PetscFree(*values);

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

205:     PetscMemcpy(*colptr,tmpcolptr,(*n+1)*sizeof(PetscInt));
206:     PetscMalloc1(((*colptr)[*n]-1),row);
207:     PetscMemcpy(*row,tmprows,((*colptr)[*n]-1)*sizeof(PetscInt));
208:     PetscMalloc1(((*colptr)[*n]-1),values);
209:     PetscMemcpy(*values,tmpvalues,((*colptr)[*n]-1)*sizeof(PetscScalar));
210:     free(tmpvalues);
211:     free(tmprows);
212:     free(tmpcolptr);

214:   }
215:   return(0);
216: }



222: /*
223:   Call clean step of PaStiX if lu->CleanUpPastix == true.
224:   Free the CSC matrix.
225:  */
226: PetscErrorCode MatDestroy_Pastix(Mat A)
227: {
228:   Mat_Pastix     *lu=(Mat_Pastix*)A->spptr;
230:   PetscMPIInt    size=lu->commSize;

233:   if (lu && lu->CleanUpPastix) {
234:     /* Terminate instance, deallocate memories */
235:     if (size > 1) {
236:       VecScatterDestroy(&lu->scat_rhs);
237:       VecDestroy(&lu->b_seq);
238:       VecScatterDestroy(&lu->scat_sol);
239:     }

241:     lu->iparm[IPARM_START_TASK]=API_TASK_CLEAN;
242:     lu->iparm[IPARM_END_TASK]  =API_TASK_CLEAN;

244:     PASTIX_CALL(&(lu->pastix_data),
245:                 lu->pastix_comm,
246:                 lu->n,
247:                 lu->colptr,
248:                 lu->row,
249:                 (PastixScalar*)lu->val,
250:                 lu->perm,
251:                 lu->invp,
252:                 (PastixScalar*)lu->rhs,
253:                 lu->rhsnbr,
254:                 lu->iparm,
255:                 lu->dparm);

257:     PetscFree(lu->colptr);
258:     PetscFree(lu->row);
259:     PetscFree(lu->val);
260:     PetscFree(lu->perm);
261:     PetscFree(lu->invp);
262:     MPI_Comm_free(&(lu->pastix_comm));
263:   }
264:   if (lu && lu->Destroy) {
265:     (lu->Destroy)(A);
266:   }
267:   PetscFree(A->spptr);
268:   return(0);
269: }

273: /*
274:   Gather right-hand-side.
275:   Call for Solve step.
276:   Scatter solution.
277:  */
278: PetscErrorCode MatSolve_PaStiX(Mat A,Vec b,Vec x)
279: {
280:   Mat_Pastix     *lu=(Mat_Pastix*)A->spptr;
281:   PetscScalar    *array;
282:   Vec            x_seq;

286:   lu->rhsnbr = 1;
287:   x_seq      = lu->b_seq;
288:   if (lu->commSize > 1) {
289:     /* PaStiX only supports centralized rhs. Scatter b into a seqential rhs vector */
290:     VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
291:     VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);
292:     VecGetArray(x_seq,&array);
293:   } else {  /* size == 1 */
294:     VecCopy(b,x);
295:     VecGetArray(x,&array);
296:   }
297:   lu->rhs = array;
298:   if (lu->commSize == 1) {
299:     VecRestoreArray(x,&array);
300:   } else {
301:     VecRestoreArray(x_seq,&array);
302:   }

304:   /* solve phase */
305:   /*-------------*/
306:   lu->iparm[IPARM_START_TASK] = API_TASK_SOLVE;
307:   lu->iparm[IPARM_END_TASK]   = API_TASK_REFINE;
308:   lu->iparm[IPARM_RHS_MAKING] = API_RHS_B;

310:   PASTIX_CALL(&(lu->pastix_data),
311:               lu->pastix_comm,
312:               lu->n,
313:               lu->colptr,
314:               lu->row,
315:               (PastixScalar*)lu->val,
316:               lu->perm,
317:               lu->invp,
318:               (PastixScalar*)lu->rhs,
319:               lu->rhsnbr,
320:               lu->iparm,
321:               lu->dparm);

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

325:   if (lu->commSize == 1) {
326:     VecRestoreArray(x,&(lu->rhs));
327:   } else {
328:     VecRestoreArray(x_seq,&(lu->rhs));
329:   }

331:   if (lu->commSize > 1) { /* convert PaStiX centralized solution to petsc mpi x */
332:     VecScatterBegin(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
333:     VecScatterEnd(lu->scat_sol,x_seq,x,INSERT_VALUES,SCATTER_FORWARD);
334:   }
335:   return(0);
336: }

338: /*
339:   Numeric factorisation using PaStiX solver.

341:  */
344: PetscErrorCode MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo *info)
345: {
346:   Mat_Pastix     *lu =(Mat_Pastix*)(F)->spptr;
347:   Mat            *tseq;
348:   PetscErrorCode 0;
349:   PetscInt       icntl;
350:   PetscInt       M=A->rmap->N;
351:   PetscBool      valOnly,flg, isSym;
352:   Mat            F_diag;
353:   IS             is_iden;
354:   Vec            b;
355:   IS             isrow;
356:   PetscBool      isSeqAIJ,isSeqSBAIJ,isMPIAIJ;

359:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
360:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
361:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
362:   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
363:     (F)->ops->solve = MatSolve_PaStiX;

365:     /* Initialize a PASTIX instance */
366:     MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->pastix_comm));
367:     MPI_Comm_rank(lu->pastix_comm, &lu->commRank);
368:     MPI_Comm_size(lu->pastix_comm, &lu->commSize);

370:     /* Set pastix options */
371:     lu->iparm[IPARM_MODIFY_PARAMETER] = API_NO;
372:     lu->iparm[IPARM_START_TASK]       = API_TASK_INIT;
373:     lu->iparm[IPARM_END_TASK]         = API_TASK_INIT;

375:     lu->rhsnbr = 1;

377:     /* Call to set default pastix options */
378:     PASTIX_CALL(&(lu->pastix_data),
379:                 lu->pastix_comm,
380:                 lu->n,
381:                 lu->colptr,
382:                 lu->row,
383:                 (PastixScalar*)lu->val,
384:                 lu->perm,
385:                 lu->invp,
386:                 (PastixScalar*)lu->rhs,
387:                 lu->rhsnbr,
388:                 lu->iparm,
389:                 lu->dparm);

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

393:     icntl = -1;

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

397:     PetscOptionsInt("-mat_pastix_verbose","iparm[IPARM_VERBOSE] : level of printing (0 to 2)","None",lu->iparm[IPARM_VERBOSE],&icntl,&flg);
398:     if ((flg && icntl >= 0) || PetscLogPrintInfo) {
399:       lu->iparm[IPARM_VERBOSE] =  icntl;
400:     }
401:     icntl=-1;
402:     PetscOptionsInt("-mat_pastix_threadnbr","iparm[IPARM_THREAD_NBR] : Number of thread by MPI node","None",lu->iparm[IPARM_THREAD_NBR],&icntl,&flg);
403:     if ((flg && icntl > 0)) {
404:       lu->iparm[IPARM_THREAD_NBR] = icntl;
405:     }
406:     PetscOptionsEnd();
407:     valOnly = PETSC_FALSE;
408:   } else {
409:     if (isSeqAIJ || isMPIAIJ) {
410:       PetscFree(lu->colptr);
411:       PetscFree(lu->row);
412:       PetscFree(lu->val);
413:       valOnly = PETSC_FALSE;
414:     } else valOnly = PETSC_TRUE;
415:   }

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

419:   /* convert mpi A to seq mat A */
420:   ISCreateStride(PETSC_COMM_SELF,M,0,1,&isrow);
421:   MatGetSubMatrices(A,1,&isrow,&isrow,MAT_INITIAL_MATRIX,&tseq);
422:   ISDestroy(&isrow);

424:   MatConvertToCSC(*tseq,valOnly, &lu->n, &lu->colptr, &lu->row, &lu->val);
425:   MatIsSymmetric(*tseq,0.0,&isSym);
426:   MatDestroyMatrices(1,&tseq);

428:   if (!lu->perm) {
429:     PetscMalloc1((lu->n),&(lu->perm));
430:     PetscMalloc1((lu->n),&(lu->invp));
431:   }

433:   if (isSym) {
434:     /* On symmetric matrix, LLT */
435:     lu->iparm[IPARM_SYM]           = API_SYM_YES;
436:     lu->iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
437:   } else {
438:     /* On unsymmetric matrix, LU */
439:     lu->iparm[IPARM_SYM]           = API_SYM_NO;
440:     lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
441:   }

443:   /*----------------*/
444:   if (lu->matstruc == DIFFERENT_NONZERO_PATTERN) {
445:     if (!(isSeqAIJ || isSeqSBAIJ)) {
446:       /* PaStiX only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */
447:       VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);
448:       ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);
449:       MatGetVecs(A,NULL,&b);
450:       VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);
451:       VecScatterCreate(lu->b_seq,is_iden,b,is_iden,&lu->scat_sol);
452:       ISDestroy(&is_iden);
453:       VecDestroy(&b);
454:     }
455:     lu->iparm[IPARM_START_TASK] = API_TASK_ORDERING;
456:     lu->iparm[IPARM_END_TASK]   = API_TASK_NUMFACT;

458:     PASTIX_CALL(&(lu->pastix_data),
459:                 lu->pastix_comm,
460:                 lu->n,
461:                 lu->colptr,
462:                 lu->row,
463:                 (PastixScalar*)lu->val,
464:                 lu->perm,
465:                 lu->invp,
466:                 (PastixScalar*)lu->rhs,
467:                 lu->rhsnbr,
468:                 lu->iparm,
469:                 lu->dparm);
470:     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]);
471:   } else {
472:     lu->iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
473:     lu->iparm[IPARM_END_TASK]   = API_TASK_NUMFACT;
474:     PASTIX_CALL(&(lu->pastix_data),
475:                 lu->pastix_comm,
476:                 lu->n,
477:                 lu->colptr,
478:                 lu->row,
479:                 (PastixScalar*)lu->val,
480:                 lu->perm,
481:                 lu->invp,
482:                 (PastixScalar*)lu->rhs,
483:                 lu->rhsnbr,
484:                 lu->iparm,
485:                 lu->dparm);

487:     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]);
488:   }

490:   if (lu->commSize > 1) {
491:     if ((F)->factortype == MAT_FACTOR_LU) {
492:       F_diag = ((Mat_MPIAIJ*)(F)->data)->A;
493:     } else {
494:       F_diag = ((Mat_MPISBAIJ*)(F)->data)->A;
495:     }
496:     F_diag->assembled = PETSC_TRUE;
497:   }
498:   (F)->assembled    = PETSC_TRUE;
499:   lu->matstruc      = SAME_NONZERO_PATTERN;
500:   lu->CleanUpPastix = PETSC_TRUE;
501:   return(0);
502: }

504: /* Note the Petsc r and c permutations are ignored */
507: PetscErrorCode MatLUFactorSymbolic_AIJPASTIX(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
508: {
509:   Mat_Pastix *lu = (Mat_Pastix*)F->spptr;

512:   lu->iparm[IPARM_FACTORIZATION] = API_FACT_LU;
513:   lu->iparm[IPARM_SYM]           = API_SYM_YES;
514:   lu->matstruc                   = DIFFERENT_NONZERO_PATTERN;
515:   F->ops->lufactornumeric        = MatFactorNumeric_PaStiX;
516:   return(0);
517: }


520: /* Note the Petsc r permutation is ignored */
523: PetscErrorCode MatCholeskyFactorSymbolic_SBAIJPASTIX(Mat F,Mat A,IS r,const MatFactorInfo *info)
524: {
525:   Mat_Pastix *lu = (Mat_Pastix*)(F)->spptr;

528:   lu->iparm[IPARM_FACTORIZATION]  = API_FACT_LLT;
529:   lu->iparm[IPARM_SYM]            = API_SYM_NO;
530:   lu->matstruc                    = DIFFERENT_NONZERO_PATTERN;
531:   (F)->ops->choleskyfactornumeric = MatFactorNumeric_PaStiX;
532:   return(0);
533: }

537: PetscErrorCode MatView_PaStiX(Mat A,PetscViewer viewer)
538: {
539:   PetscErrorCode    ierr;
540:   PetscBool         iascii;
541:   PetscViewerFormat format;

544:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
545:   if (iascii) {
546:     PetscViewerGetFormat(viewer,&format);
547:     if (format == PETSC_VIEWER_ASCII_INFO) {
548:       Mat_Pastix *lu=(Mat_Pastix*)A->spptr;

550:       PetscViewerASCIIPrintf(viewer,"PaStiX run parameters:\n");
551:       PetscViewerASCIIPrintf(viewer,"  Matrix type :                      %s \n",((lu->iparm[IPARM_SYM] == API_SYM_YES) ? "Symmetric" : "Unsymmetric"));
552:       PetscViewerASCIIPrintf(viewer,"  Level of printing (0,1,2):         %d \n",lu->iparm[IPARM_VERBOSE]);
553:       PetscViewerASCIIPrintf(viewer,"  Number of refinements iterations : %d \n",lu->iparm[IPARM_NBITER]);
554:       PetscPrintf(PETSC_COMM_SELF,"  Error :                        %g \n",lu->dparm[DPARM_RELATIVE_ERROR]);
555:     }
556:   }
557:   return(0);
558: }


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

565:   Use ./configure --download-pastix to have PETSc installed with PaStiX

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

571:   Level: beginner

573: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage

575: M*/


580: PetscErrorCode MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo *info)
581: {
582:   Mat_Pastix *lu =(Mat_Pastix*)A->spptr;

585:   info->block_size        = 1.0;
586:   info->nz_allocated      = lu->iparm[IPARM_NNZEROS];
587:   info->nz_used           = lu->iparm[IPARM_NNZEROS];
588:   info->nz_unneeded       = 0.0;
589:   info->assemblies        = 0.0;
590:   info->mallocs           = 0.0;
591:   info->memory            = 0.0;
592:   info->fill_ratio_given  = 0;
593:   info->fill_ratio_needed = 0;
594:   info->factor_mallocs    = 0;
595:   return(0);
596: }

600: PetscErrorCode MatFactorGetSolverPackage_pastix(Mat A,const MatSolverPackage *type)
601: {
603:   *type = MATSOLVERPASTIX;
604:   return(0);
605: }

607: /*
608:     The seq and mpi versions of this function are the same
609: */
612: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat *F)
613: {
614:   Mat            B;
616:   Mat_Pastix     *pastix;

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

626:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
627:   B->ops->view             = MatView_PaStiX;
628:   B->ops->getinfo          = MatGetInfo_PaStiX;

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

632:   B->factortype = MAT_FACTOR_LU;

634:   PetscNewLog(B,&pastix);

636:   pastix->CleanUpPastix = PETSC_FALSE;
637:   pastix->isAIJ         = PETSC_TRUE;
638:   pastix->scat_rhs      = NULL;
639:   pastix->scat_sol      = NULL;
640:   pastix->Destroy       = B->ops->destroy;
641:   B->ops->destroy       = MatDestroy_Pastix;
642:   B->spptr              = (void*)pastix;

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

650: PETSC_EXTERN PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat *F)
651: {
652:   Mat            B;
654:   Mat_Pastix     *pastix;

657:   if (ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
658:   /* Create the factorization matrix */
659:   MatCreate(PetscObjectComm((PetscObject)A),&B);
660:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
661:   MatSetType(B,((PetscObject)A)->type_name);
662:   MatSeqAIJSetPreallocation(B,0,NULL);
663:   MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);

665:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJPASTIX;
666:   B->ops->view             = MatView_PaStiX;

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

670:   B->factortype = MAT_FACTOR_LU;

672:   PetscNewLog(B,&pastix);

674:   pastix->CleanUpPastix = PETSC_FALSE;
675:   pastix->isAIJ         = PETSC_TRUE;
676:   pastix->scat_rhs      = NULL;
677:   pastix->scat_sol      = NULL;
678:   pastix->Destroy       = B->ops->destroy;
679:   B->ops->destroy       = MatDestroy_Pastix;
680:   B->spptr              = (void*)pastix;

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

688: PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
689: {
690:   Mat            B;
692:   Mat_Pastix     *pastix;

695:   if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
696:   /* Create the factorization matrix */
697:   MatCreate(PetscObjectComm((PetscObject)A),&B);
698:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
699:   MatSetType(B,((PetscObject)A)->type_name);
700:   MatSeqSBAIJSetPreallocation(B,1,0,NULL);
701:   MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);

703:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
704:   B->ops->view                   = MatView_PaStiX;

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

708:   B->factortype = MAT_FACTOR_CHOLESKY;

710:   PetscNewLog(B,&pastix);

712:   pastix->CleanUpPastix = PETSC_FALSE;
713:   pastix->isAIJ         = PETSC_TRUE;
714:   pastix->scat_rhs      = NULL;
715:   pastix->scat_sol      = NULL;
716:   pastix->Destroy       = B->ops->destroy;
717:   B->ops->destroy       = MatDestroy_Pastix;
718:   B->spptr              = (void*)pastix;

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

726: PETSC_EXTERN PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
727: {
728:   Mat            B;
730:   Mat_Pastix     *pastix;

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

735:   /* Create the factorization matrix */
736:   MatCreate(PetscObjectComm((PetscObject)A),&B);
737:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
738:   MatSetType(B,((PetscObject)A)->type_name);
739:   MatSeqSBAIJSetPreallocation(B,1,0,NULL);
740:   MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);

742:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
743:   B->ops->view                   = MatView_PaStiX;

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

747:   B->factortype = MAT_FACTOR_CHOLESKY;

749:   PetscNewLog(B,&pastix);

751:   pastix->CleanUpPastix = PETSC_FALSE;
752:   pastix->isAIJ         = PETSC_TRUE;
753:   pastix->scat_rhs      = NULL;
754:   pastix->scat_sol      = NULL;
755:   pastix->Destroy       = B->ops->destroy;
756:   B->ops->destroy       = MatDestroy_Pastix;
757:   B->spptr              = (void*)pastix;

759:   *F = B;
760:   return(0);
761: }