Actual source code: superlu.c

      1: #define PETSCMAT_DLL

      3: /*
      4:         Provides an interface to the SuperLU 3.0 sparse solver
      5: */

     7:  #include src/mat/impls/aij/seq/aij.h

     10: #if defined(PETSC_USE_COMPLEX)
     11: #include "zsp_defs.h"
     12: #else
     13: #include "dsp_defs.h"
     14: #endif
     15: #include "util.h"

     18: typedef struct {
     19:   SuperMatrix       A,L,U,B,X;
     20:   superlu_options_t options;
     21:   PetscInt          *perm_c; /* column permutation vector */
     22:   PetscInt          *perm_r; /* row permutations from partial pivoting */
     23:   PetscInt          *etree;
     24:   PetscReal         *R, *C;
     25:   char              equed[1];
     26:   PetscInt          lwork;
     27:   void              *work;
     28:   PetscReal         rpg, rcond;
     29:   mem_usage_t       mem_usage;
     30:   MatStructure      flg;

     32:   /* A few function pointers for inheritance */
     33:   PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
     34:   PetscErrorCode (*MatView)(Mat,PetscViewer);
     35:   PetscErrorCode (*MatAssemblyEnd)(Mat,MatAssemblyType);
     36:   PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
     37:   PetscErrorCode (*MatDestroy)(Mat);

     39:   /* Flag to clean up (non-global) SuperLU objects during Destroy */
     40:   PetscTruth CleanUpSuperLU;
     41: } Mat_SuperLU;

     44: EXTERN PetscErrorCode MatFactorInfo_SuperLU(Mat,PetscViewer);
     45: EXTERN PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat,IS,IS,MatFactorInfo*,Mat*);

     48: EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SuperLU_SeqAIJ(Mat,MatType,MatReuse,Mat*);
     49: EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SuperLU(Mat,MatType,MatReuse,Mat*);

     54: PetscErrorCode MatDestroy_SuperLU(Mat A)
     55: {
     57:   Mat_SuperLU    *lu=(Mat_SuperLU*)A->spptr;

     60:   if (lu->CleanUpSuperLU) { /* Free the SuperLU datastructures */
     61:     Destroy_SuperMatrix_Store(&lu->A);
     62:     Destroy_SuperMatrix_Store(&lu->B);
     63:     Destroy_SuperMatrix_Store(&lu->X);

     65:     PetscFree(lu->etree);
     66:     PetscFree(lu->perm_r);
     67:     PetscFree(lu->perm_c);
     68:     PetscFree(lu->R);
     69:     PetscFree(lu->C);
     70:     if ( lu->lwork >= 0 ) {
     71:       Destroy_SuperNode_Matrix(&lu->L);
     72:       Destroy_CompCol_Matrix(&lu->U);
     73:     }
     74:   }
     75:   MatConvert_SuperLU_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);
     76:   (*A->ops->destroy)(A);
     77:   return(0);
     78: }

     82: PetscErrorCode MatView_SuperLU(Mat A,PetscViewer viewer)
     83: {
     84:   PetscErrorCode    ierr;
     85:   PetscTruth        iascii;
     86:   PetscViewerFormat format;
     87:   Mat_SuperLU       *lu=(Mat_SuperLU*)(A->spptr);

     90:   (*lu->MatView)(A,viewer);

     92:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     93:   if (iascii) {
     94:     PetscViewerGetFormat(viewer,&format);
     95:     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
     96:       MatFactorInfo_SuperLU(A,viewer);
     97:     }
     98:   }
     99:   return(0);
    100: }

    104: PetscErrorCode MatAssemblyEnd_SuperLU(Mat A,MatAssemblyType mode) {
    106:   Mat_SuperLU    *lu=(Mat_SuperLU*)(A->spptr);

    109:   (*lu->MatAssemblyEnd)(A,mode);
    110:   lu->MatLUFactorSymbolic  = A->ops->lufactorsymbolic;
    111:   A->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU;
    112:   return(0);
    113: }

    115: /* This function was written for SuperLU 2.0 by Matthew Knepley. Not tested for SuperLU 3.0! */
    116: #ifdef SuperLU2
    117:  #include src/mat/impls/dense/seq/dense.h
    120: PetscErrorCode MatCreateNull_SuperLU(Mat A,Mat *nullMat)
    121: {
    122:   Mat_SuperLU   *lu = (Mat_SuperLU*)A->spptr;
    123:   PetscInt      numRows = A->m,numCols = A->n;
    124:   SCformat      *Lstore;
    125:   PetscInt      numNullCols,size;
    126:   SuperLUStat_t stat;
    127: #if defined(PETSC_USE_COMPLEX)
    128:   doublecomplex *nullVals,*workVals;
    129: #else
    130:   PetscScalar   *nullVals,*workVals;
    131: #endif
    132:   PetscInt           row,newRow,col,newCol,block,b;

    136:   if (!A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    137:   numNullCols = numCols - numRows;
    138:   if (numNullCols < 0) SETERRQ(PETSC_ERR_ARG_WRONG,"Function only applies to underdetermined problems");
    139:   /* Create the null matrix using MATSEQDENSE explicitly */
    140:   MatCreate(A->comm,nullMat);
    141:   MatSetSizes(*nullMat,numRows,numNullCols,numRows,numNullCols);
    142:   MatSetType(*nullMat,MATSEQDENSE);
    143:   MatSeqDenseSetPreallocation(*nullMat,PETSC_NULL);
    144:   if (!numNullCols) {
    145:     MatAssemblyBegin(*nullMat,MAT_FINAL_ASSEMBLY);
    146:     MatAssemblyEnd(*nullMat,MAT_FINAL_ASSEMBLY);
    147:     return(0);
    148:   }
    149: #if defined(PETSC_USE_COMPLEX)
    150:   nullVals = (doublecomplex*)((Mat_SeqDense*)(*nullMat)->data)->v;
    151: #else
    152:   nullVals = ((Mat_SeqDense*)(*nullMat)->data)->v;
    153: #endif
    154:   /* Copy in the columns */
    155:   Lstore = (SCformat*)lu->L.Store;
    156:   for(block = 0; block <= Lstore->nsuper; block++) {
    157:     newRow = Lstore->sup_to_col[block];
    158:     size   = Lstore->sup_to_col[block+1] - Lstore->sup_to_col[block];
    159:     for(col = Lstore->rowind_colptr[newRow]; col < Lstore->rowind_colptr[newRow+1]; col++) {
    160:       newCol = Lstore->rowind[col];
    161:       if (newCol >= numRows) {
    162:         for(b = 0; b < size; b++)
    163: #if defined(PETSC_USE_COMPLEX)
    164:           nullVals[(newCol-numRows)*numRows+newRow+b] = ((doublecomplex*)Lstore->nzval)[Lstore->nzval_colptr[newRow+b]+col];
    165: #else
    166:           nullVals[(newCol-numRows)*numRows+newRow+b] = ((double*)Lstore->nzval)[Lstore->nzval_colptr[newRow+b]+col];
    167: #endif
    168:       }
    169:     }
    170:   }
    171:   /* Permute rhs to form P^T_c B */
    172:   PetscMalloc(numRows*sizeof(PetscReal),&workVals);
    173:   for(b = 0; b < numNullCols; b++) {
    174:     for(row = 0; row < numRows; row++) workVals[lu->perm_c[row]] = nullVals[b*numRows+row];
    175:     for(row = 0; row < numRows; row++) nullVals[b*numRows+row]   = workVals[row];
    176:   }
    177:   /* Backward solve the upper triangle A x = b */
    178:   for(b = 0; b < numNullCols; b++) {
    179: #if defined(PETSC_USE_COMPLEX)
    180:     sp_ztrsv("L","T","U",&lu->L,&lu->U,&nullVals[b*numRows],&stat,&ierr);
    181: #else
    182:     sp_dtrsv("L","T","U",&lu->L,&lu->U,&nullVals[b*numRows],&stat,&ierr);
    183: #endif
    184:     if (ierr < 0)
    185:       SETERRQ1(PETSC_ERR_ARG_WRONG,"The argument %D was invalid",-ierr);
    186:   }
    187:   PetscFree(workVals);

    189:   MatAssemblyBegin(*nullMat,MAT_FINAL_ASSEMBLY);
    190:   MatAssemblyEnd(*nullMat,MAT_FINAL_ASSEMBLY);
    191:   return(0);
    192: }
    193: #endif

    197: PetscErrorCode MatSolve_SuperLU(Mat A,Vec b,Vec x)
    198: {
    199:   Mat_SuperLU    *lu = (Mat_SuperLU*)A->spptr;
    200:   PetscScalar    *barray,*xarray;
    202:   PetscInt       info,i;
    203:   SuperLUStat_t  stat;
    204:   PetscReal      ferr,berr;

    207:   if ( lu->lwork == -1 ) {
    208:     return(0);
    209:   }
    210:   lu->B.ncol = 1;   /* Set the number of right-hand side */
    211:   VecGetArray(b,&barray);
    212:   VecGetArray(x,&xarray);

    214: #if defined(PETSC_USE_COMPLEX)
    215:   ((DNformat*)lu->B.Store)->nzval = (doublecomplex*)barray;
    216:   ((DNformat*)lu->X.Store)->nzval = (doublecomplex*)xarray;
    217: #else
    218:   ((DNformat*)lu->B.Store)->nzval = barray;
    219:   ((DNformat*)lu->X.Store)->nzval = xarray;
    220: #endif

    222:   /* Initialize the statistics variables. */
    223:   StatInit(&stat);

    225:   lu->options.Fact  = FACTORED; /* Indicate the factored form of A is supplied. */
    226:   lu->options.Trans = TRANS;
    227: #if defined(PETSC_USE_COMPLEX)
    228:   zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
    229:            &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
    230:            &lu->mem_usage, &stat, &info);
    231: #else
    232:   dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
    233:            &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
    234:            &lu->mem_usage, &stat, &info);
    235: #endif
    236:   VecRestoreArray(b,&barray);
    237:   VecRestoreArray(x,&xarray);

    239:   if ( !info || info == lu->A.ncol+1 ) {
    240:     if ( lu->options.IterRefine ) {
    241:       PetscPrintf(PETSC_COMM_SELF,"Iterative Refinement:\n");
    242:       PetscPrintf(PETSC_COMM_SELF,"  %8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR");
    243:       for (i = 0; i < 1; ++i)
    244:         PetscPrintf(PETSC_COMM_SELF,"  %8d%8d%16e%16e\n", i+1, stat.RefineSteps, ferr, berr);
    245:     }
    246:   } else if ( info > 0 ){
    247:     if ( lu->lwork == -1 ) {
    248:       PetscPrintf(PETSC_COMM_SELF,"  ** Estimated memory: %D bytes\n", info - lu->A.ncol);
    249:     } else {
    250:       PetscPrintf(PETSC_COMM_SELF,"  Warning: gssvx() returns info %D\n",info);
    251:     }
    252:   } else if (info < 0){
    253:     SETERRQ2(PETSC_ERR_LIB, "info = %D, the %D-th argument in gssvx() had an illegal value", info,-info);
    254:   }

    256:   if ( lu->options.PrintStat ) {
    257:     PetscPrintf(PETSC_COMM_SELF,"MatSolve__SuperLU():\n");
    258:     StatPrint(&stat);
    259:   }
    260:   StatFree(&stat);
    261:   return(0);
    262: }

    266: PetscErrorCode MatLUFactorNumeric_SuperLU(Mat A,MatFactorInfo *info,Mat *F)
    267: {
    268:   Mat_SeqAIJ     *aa = (Mat_SeqAIJ*)(A)->data;
    269:   Mat_SuperLU    *lu = (Mat_SuperLU*)(*F)->spptr;
    271:   PetscInt       sinfo;
    272:   SuperLUStat_t  stat;
    273:   PetscReal      ferr, berr;
    274:   NCformat       *Ustore;
    275:   SCformat       *Lstore;
    276:
    278:   if (lu->flg == SAME_NONZERO_PATTERN){ /* successing numerical factorization */
    279:     lu->options.Fact = SamePattern;
    280:     /* Ref: ~SuperLU_3.0/EXAMPLE/dlinsolx2.c */
    281:     Destroy_SuperMatrix_Store(&lu->A);
    282:     if ( lu->lwork >= 0 ) {
    283:       Destroy_SuperNode_Matrix(&lu->L);
    284:       Destroy_CompCol_Matrix(&lu->U);
    285:       lu->options.Fact = SamePattern;
    286:     }
    287:   }

    289:   /* Create the SuperMatrix for lu->A=A^T:
    290:        Since SuperLU likes column-oriented matrices,we pass it the transpose,
    291:        and then solve A^T X = B in MatSolve(). */
    292: #if defined(PETSC_USE_COMPLEX)
    293:   zCreate_CompCol_Matrix(&lu->A,A->n,A->m,aa->nz,(doublecomplex*)aa->a,aa->j,aa->i,
    294:                            SLU_NC,SLU_Z,SLU_GE);
    295: #else
    296:   dCreate_CompCol_Matrix(&lu->A,A->n,A->m,aa->nz,aa->a,aa->j,aa->i,
    297:                            SLU_NC,SLU_D,SLU_GE);
    298: #endif
    299:
    300:   /* Initialize the statistics variables. */
    301:   StatInit(&stat);

    303:   /* Numerical factorization */
    304:   lu->B.ncol = 0;  /* Indicate not to solve the system */
    305: #if defined(PETSC_USE_COMPLEX)
    306:    zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
    307:            &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
    308:            &lu->mem_usage, &stat, &sinfo);
    309: #else
    310:   dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
    311:            &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
    312:            &lu->mem_usage, &stat, &sinfo);
    313: #endif
    314:   if ( !sinfo || sinfo == lu->A.ncol+1 ) {
    315:     if ( lu->options.PivotGrowth )
    316:       PetscPrintf(PETSC_COMM_SELF,"  Recip. pivot growth = %e\n", lu->rpg);
    317:     if ( lu->options.ConditionNumber )
    318:       PetscPrintf(PETSC_COMM_SELF,"  Recip. condition number = %e\n", lu->rcond);
    319:   } else if ( sinfo > 0 ){
    320:     if ( lu->lwork == -1 ) {
    321:       PetscPrintf(PETSC_COMM_SELF,"  ** Estimated memory: %D bytes\n", sinfo - lu->A.ncol);
    322:     } else {
    323:       PetscPrintf(PETSC_COMM_SELF,"  Warning: gssvx() returns info %D\n",sinfo);
    324:     }
    325:   } else { /* sinfo < 0 */
    326:     SETERRQ2(PETSC_ERR_LIB, "info = %D, the %D-th argument in gssvx() had an illegal value", sinfo,-sinfo);
    327:   }

    329:   if ( lu->options.PrintStat ) {
    330:     PetscPrintf(PETSC_COMM_SELF,"MatLUFactorNumeric_SuperLU():\n");
    331:     StatPrint(&stat);
    332:     Lstore = (SCformat *) lu->L.Store;
    333:     Ustore = (NCformat *) lu->U.Store;
    334:     PetscPrintf(PETSC_COMM_SELF,"  No of nonzeros in factor L = %D\n", Lstore->nnz);
    335:     PetscPrintf(PETSC_COMM_SELF,"  No of nonzeros in factor U = %D\n", Ustore->nnz);
    336:     PetscPrintf(PETSC_COMM_SELF,"  No of nonzeros in L+U = %D\n", Lstore->nnz + Ustore->nnz - lu->A.ncol);
    337:     PetscPrintf(PETSC_COMM_SELF,"  L\\U MB %.3f\ttotal MB needed %.3f\texpansions %D\n",
    338:                lu->mem_usage.for_lu/1e6, lu->mem_usage.total_needed/1e6,
    339:                lu->mem_usage.expansions);
    340:   }
    341:   StatFree(&stat);

    343:   lu->flg = SAME_NONZERO_PATTERN;
    344:   return(0);
    345: }

    347: /*
    348:    Note the r permutation is ignored
    349: */
    352: PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
    353: {
    354:   Mat            B;
    355:   Mat_SuperLU    *lu;
    357:   PetscInt       m=A->m,n=A->n,indx;
    358:   PetscTruth     flg;
    359:   const char   *colperm[]={"NATURAL","MMD_ATA","MMD_AT_PLUS_A","COLAMD"}; /* MY_PERMC - not supported by the petsc interface yet */
    360:   const char   *iterrefine[]={"NOREFINE", "SINGLE", "DOUBLE", "EXTRA"};
    361:   const char   *rowperm[]={"NOROWPERM", "LargeDiag"}; /* MY_PERMC - not supported by the petsc interface yet */

    364:   MatCreate(A->comm,&B);
    365:   MatSetSizes(B,A->m,A->n,PETSC_DETERMINE,PETSC_DETERMINE);
    366:   MatSetType(B,A->type_name);
    367:   MatSeqAIJSetPreallocation(B,0,PETSC_NULL);

    369:   B->ops->lufactornumeric = MatLUFactorNumeric_SuperLU;
    370:   B->ops->solve           = MatSolve_SuperLU;
    371:   B->factor               = FACTOR_LU;
    372:   B->assembled            = PETSC_TRUE;  /* required by -ksp_view */
    373:
    374:   lu = (Mat_SuperLU*)(B->spptr);

    376:   /* Set SuperLU options */
    377:     /* the default values for options argument:
    378:         options.Fact = DOFACT;
    379:         options.Equil = YES;
    380:             options.ColPerm = COLAMD;
    381:         options.DiagPivotThresh = 1.0;
    382:             options.Trans = NOTRANS;
    383:             options.IterRefine = NOREFINE;
    384:             options.SymmetricMode = NO;
    385:             options.PivotGrowth = NO;
    386:             options.ConditionNumber = NO;
    387:             options.PrintStat = YES;
    388:     */
    389:   set_default_options(&lu->options);
    390:   /* equilibration causes error in solve(), thus not supported here. See dgssvx.c for possible reason. */
    391:   lu->options.Equil = NO;
    392:   lu->options.PrintStat = NO;
    393:   lu->lwork = 0;   /* allocate space internally by system malloc */

    395:   PetscOptionsBegin(A->comm,A->prefix,"SuperLU Options","Mat");
    396:   /*
    397:   PetscOptionsTruth("-mat_superlu_equil","Equil","None",PETSC_FALSE,&flg,0);
    398:   if (flg) lu->options.Equil = YES; -- not supported by the interface !!!
    399:   */
    400:   PetscOptionsEList("-mat_superlu_colperm","ColPerm","None",colperm,4,colperm[3],&indx,&flg);
    401:   if (flg) {lu->options.ColPerm = (colperm_t)indx;}
    402:   PetscOptionsEList("-mat_superlu_iterrefine","IterRefine","None",iterrefine,4,iterrefine[0],&indx,&flg);
    403:   if (flg) { lu->options.IterRefine = (IterRefine_t)indx;}
    404:   PetscOptionsTruth("-mat_superlu_symmetricmode","SymmetricMode","None",PETSC_FALSE,&flg,0);
    405:   if (flg) lu->options.SymmetricMode = YES;
    406:   PetscOptionsReal("-mat_superlu_diagpivotthresh","DiagPivotThresh","None",lu->options.DiagPivotThresh,&lu->options.DiagPivotThresh,PETSC_NULL);
    407:   PetscOptionsTruth("-mat_superlu_pivotgrowth","PivotGrowth","None",PETSC_FALSE,&flg,0);
    408:   if (flg) lu->options.PivotGrowth = YES;
    409:   PetscOptionsTruth("-mat_superlu_conditionnumber","ConditionNumber","None",PETSC_FALSE,&flg,0);
    410:   if (flg) lu->options.ConditionNumber = YES;
    411:   PetscOptionsEList("-mat_superlu_rowperm","rowperm","None",rowperm,2,rowperm[0],&indx,&flg);
    412:   if (flg) {lu->options.RowPerm = (rowperm_t)indx;}
    413:   PetscOptionsTruth("-mat_superlu_replacetinypivot","ReplaceTinyPivot","None",PETSC_FALSE,&flg,0);
    414:   if (flg) lu->options.ReplaceTinyPivot = YES;
    415:   PetscOptionsTruth("-mat_superlu_printstat","PrintStat","None",PETSC_FALSE,&flg,0);
    416:   if (flg) lu->options.PrintStat = YES;
    417:   PetscOptionsInt("-mat_superlu_lwork","size of work array in bytes used by factorization","None",lu->lwork,&lu->lwork,PETSC_NULL);
    418:   if (lu->lwork > 0 ){
    419:     PetscMalloc(lu->lwork,&lu->work);
    420:   } else if (lu->lwork != 0 && lu->lwork != -1){
    421:     PetscPrintf(PETSC_COMM_SELF,"   Warning: lwork %D is not supported by SUPERLU. The default lwork=0 is used.\n",lu->lwork);
    422:     lu->lwork = 0;
    423:   }
    424:   PetscOptionsEnd();

    426: #ifdef SUPERLU2
    427:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatCreateNull","MatCreateNull_SuperLU",
    428:                                     (void(*)(void))MatCreateNull_SuperLU);
    429: #endif

    431:   /* Allocate spaces (notice sizes are for the transpose) */
    432:   PetscMalloc(m*sizeof(PetscInt),&lu->etree);
    433:   PetscMalloc(n*sizeof(PetscInt),&lu->perm_r);
    434:   PetscMalloc(m*sizeof(PetscInt),&lu->perm_c);
    435:   PetscMalloc(n*sizeof(PetscInt),&lu->R);
    436:   PetscMalloc(m*sizeof(PetscInt),&lu->C);
    437:
    438:   /* create rhs and solution x without allocate space for .Store */
    439: #if defined(PETSC_USE_COMPLEX)
    440:   zCreate_Dense_Matrix(&lu->B, m, 1, PETSC_NULL, m, SLU_DN, SLU_Z, SLU_GE);
    441:   zCreate_Dense_Matrix(&lu->X, m, 1, PETSC_NULL, m, SLU_DN, SLU_Z, SLU_GE);
    442: #else
    443:   dCreate_Dense_Matrix(&lu->B, m, 1, PETSC_NULL, m, SLU_DN, SLU_D, SLU_GE);
    444:   dCreate_Dense_Matrix(&lu->X, m, 1, PETSC_NULL, m, SLU_DN, SLU_D, SLU_GE);
    445: #endif

    447:   lu->flg            = DIFFERENT_NONZERO_PATTERN;
    448:   lu->CleanUpSuperLU = PETSC_TRUE;

    450:   *F = B;
    451:   PetscLogObjectMemory(B,(A->m+A->n)*sizeof(PetscInt)+sizeof(Mat_SuperLU));
    452:   return(0);
    453: }

    455: /* used by -ksp_view */
    458: PetscErrorCode MatFactorInfo_SuperLU(Mat A,PetscViewer viewer)
    459: {
    460:   Mat_SuperLU       *lu= (Mat_SuperLU*)A->spptr;
    461:   PetscErrorCode    ierr;
    462:   superlu_options_t options;

    465:   /* check if matrix is superlu_dist type */
    466:   if (A->ops->solve != MatSolve_SuperLU) return(0);

    468:   options = lu->options;
    469:   PetscViewerASCIIPrintf(viewer,"SuperLU run parameters:\n");
    470:   PetscViewerASCIIPrintf(viewer,"  Equil: %s\n",(options.Equil != NO) ? "YES": "NO");
    471:   PetscViewerASCIIPrintf(viewer,"  ColPerm: %D\n",options.ColPerm);
    472:   PetscViewerASCIIPrintf(viewer,"  IterRefine: %D\n",options.IterRefine);
    473:   PetscViewerASCIIPrintf(viewer,"  SymmetricMode: %s\n",(options.SymmetricMode != NO) ? "YES": "NO");
    474:   PetscViewerASCIIPrintf(viewer,"  DiagPivotThresh: %g\n",options.DiagPivotThresh);
    475:   PetscViewerASCIIPrintf(viewer,"  PivotGrowth: %s\n",(options.PivotGrowth != NO) ? "YES": "NO");
    476:   PetscViewerASCIIPrintf(viewer,"  ConditionNumber: %s\n",(options.ConditionNumber != NO) ? "YES": "NO");
    477:   PetscViewerASCIIPrintf(viewer,"  RowPerm: %D\n",options.RowPerm);
    478:   PetscViewerASCIIPrintf(viewer,"  ReplaceTinyPivot: %s\n",(options.ReplaceTinyPivot != NO) ? "YES": "NO");
    479:   PetscViewerASCIIPrintf(viewer,"  PrintStat: %s\n",(options.PrintStat != NO) ? "YES": "NO");
    480:   PetscViewerASCIIPrintf(viewer,"  lwork: %D\n",lu->lwork);

    482:   return(0);
    483: }

    487: PetscErrorCode MatDuplicate_SuperLU(Mat A, MatDuplicateOption op, Mat *M) {
    489:   Mat_SuperLU    *lu=(Mat_SuperLU *)A->spptr;

    492:   (*lu->MatDuplicate)(A,op,M);
    493:   PetscMemcpy((*M)->spptr,lu,sizeof(Mat_SuperLU));
    494:   return(0);
    495: }

    500: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SuperLU_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat)
    501: {
    502:   /* This routine is only called to convert an unfactored PETSc-SuperLU matrix */
    503:   /* to its base PETSc type, so we will ignore 'MatType type'. */
    505:   Mat            B=*newmat;
    506:   Mat_SuperLU    *lu=(Mat_SuperLU *)A->spptr;

    509:   if (reuse == MAT_INITIAL_MATRIX) {
    510:     MatDuplicate(A,MAT_COPY_VALUES,&B);
    511:   }
    512:   /* Reset the original function pointers */
    513:   B->ops->duplicate        = lu->MatDuplicate;
    514:   B->ops->view             = lu->MatView;
    515:   B->ops->assemblyend      = lu->MatAssemblyEnd;
    516:   B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic;
    517:   B->ops->destroy          = lu->MatDestroy;
    518:   /* lu is only a function pointer stash unless we've factored the matrix, which we haven't! */
    519:   PetscFree(lu);

    521:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_superlu_C","",PETSC_NULL);
    522:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_superlu_seqaij_C","",PETSC_NULL);

    524:   PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
    525:   *newmat = B;
    526:   return(0);
    527: }

    533: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SuperLU(Mat A,MatType type,MatReuse reuse,Mat *newmat)
    534: {
    535:   /* This routine is only called to convert to MATSUPERLU */
    536:   /* from MATSEQAIJ, so we will ignore 'MatType type'. */
    538:   Mat            B=*newmat;
    539:   Mat_SuperLU    *lu;

    542:   if (reuse == MAT_INITIAL_MATRIX) {
    543:     MatDuplicate(A,MAT_COPY_VALUES,&B);
    544:   }

    546:   PetscNew(Mat_SuperLU,&lu);
    547:   lu->MatDuplicate         = A->ops->duplicate;
    548:   lu->MatView              = A->ops->view;
    549:   lu->MatAssemblyEnd       = A->ops->assemblyend;
    550:   lu->MatLUFactorSymbolic  = A->ops->lufactorsymbolic;
    551:   lu->MatDestroy           = A->ops->destroy;
    552:   lu->CleanUpSuperLU       = PETSC_FALSE;

    554:   B->spptr                 = (void*)lu;
    555:   B->ops->duplicate        = MatDuplicate_SuperLU;
    556:   B->ops->view             = MatView_SuperLU;
    557:   B->ops->assemblyend      = MatAssemblyEnd_SuperLU;
    558:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU;
    559:   B->ops->choleskyfactorsymbolic = 0;
    560:   B->ops->destroy          = MatDestroy_SuperLU;

    562:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_superlu_C",
    563:                                            "MatConvert_SeqAIJ_SuperLU",MatConvert_SeqAIJ_SuperLU);
    564:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_superlu_seqaij_C",
    565:                                            "MatConvert_SuperLU_SeqAIJ",MatConvert_SuperLU_SeqAIJ);
    566:   PetscLogInfo((0,"MatConvert_SeqAIJ_SuperLU:Using SuperLU for SeqAIJ LU factorization and solves.\n"));
    567:   PetscObjectChangeTypeName((PetscObject)B,MATSUPERLU);
    568:   *newmat = B;
    569:   return(0);
    570: }

    573: /*MC
    574:   MATSUPERLU - MATSUPERLU = "superlu" - A matrix type providing direct solvers (LU) for sequential matrices
    575:   via the external package SuperLU.

    577:   If SuperLU is installed (see the manual for
    578:   instructions on how to declare the existence of external packages),
    579:   a matrix type can be constructed which invokes SuperLU solvers.
    580:   After calling MatCreate(...,A), simply call MatSetType(A,MATSUPERLU).

    582:   This matrix inherits from MATSEQAIJ.  As a result, MatSeqAIJSetPreallocation is
    583:   supported for this matrix type.  One can also call MatConvert for an inplace conversion to or from
    584:   the MATSEQAIJ type without data copy.

    586:   Options Database Keys:
    587: + -mat_type superlu - sets the matrix type to "superlu" during a call to MatSetFromOptions()
    588: - -mat_superlu_ordering <0,1,2,3> - 0: natural ordering,
    589:                                     1: MMD applied to A'*A,
    590:                                     2: MMD applied to A'+A,
    591:                                     3: COLAMD, approximate minimum degree column ordering

    593:    Level: beginner

    595: .seealso: PCLU
    596: M*/

    601: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SuperLU(Mat A)
    602: {

    606:   /* Change type name before calling MatSetType to force proper construction of SeqAIJ and SUPERLU types */
    607:   PetscObjectChangeTypeName((PetscObject)A,MATSUPERLU);
    608:   MatSetType(A,MATSEQAIJ);
    609:   MatConvert_SeqAIJ_SuperLU(A,MATSUPERLU,MAT_REUSE_MATRIX,&A);
    610:   return(0);
    611: }