Actual source code: factimpl.c

      2: #include <../src/ksp/pc/impls/factor/factor.h>

      4: /* ------------------------------------------------------------------------------------------*/

      7: PetscErrorCode PCFactorSetUpMatSolverType_Factor(PC pc)
      8: {
      9:   PC_Factor      *icc = (PC_Factor*)pc->data;

     13:   if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"You can only call this routine after the matrix object has been provided to the solver, for example with KSPSetOperators() or SNESSetJacobian()");
     14:   if (!pc->setupcalled && !((PC_Factor*)icc)->fact) {
     15:     MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,((PC_Factor*)icc)->factortype,&((PC_Factor*)icc)->fact);
     16:   }
     17:   return(0);
     18: }

     20: PetscErrorCode  PCFactorSetZeroPivot_Factor(PC pc,PetscReal z)
     21: {
     22:   PC_Factor *ilu = (PC_Factor*)pc->data;

     25:   ilu->info.zeropivot = z;
     26:   return(0);
     27: }

     29: PetscErrorCode  PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype)
     30: {
     31:   PC_Factor *dir = (PC_Factor*)pc->data;

     34:   if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
     35:   else {
     36:     dir->info.shifttype = (PetscReal) shifttype;
     37:     if ((shifttype == MAT_SHIFT_NONZERO || shifttype ==  MAT_SHIFT_INBLOCKS) && dir->info.shiftamount == 0.0) {
     38:       dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON; /* set default amount if user has not called PCFactorSetShiftAmount() yet */
     39:     }
     40:   }
     41:   return(0);
     42: }

     44: PetscErrorCode  PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount)
     45: {
     46:   PC_Factor *dir = (PC_Factor*)pc->data;

     49:   if (shiftamount == (PetscReal) PETSC_DECIDE) dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
     50:   else dir->info.shiftamount = shiftamount;
     51:   return(0);
     52: }

     54: PetscErrorCode  PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
     55: {
     56:   PC_Factor *ilu = (PC_Factor*)pc->data;

     59:   if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
     60:     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change tolerance after use");
     61:   }
     62:   ilu->info.usedt   = PETSC_TRUE;
     63:   ilu->info.dt      = dt;
     64:   ilu->info.dtcol   = dtcol;
     65:   ilu->info.dtcount = dtcount;
     66:   ilu->info.fill    = PETSC_DEFAULT;
     67:   return(0);
     68: }

     70: PetscErrorCode  PCFactorSetFill_Factor(PC pc,PetscReal fill)
     71: {
     72:   PC_Factor *dir = (PC_Factor*)pc->data;

     75:   dir->info.fill = fill;
     76:   return(0);
     77: }

     79: PetscErrorCode  PCFactorSetMatOrderingType_Factor(PC pc,MatOrderingType ordering)
     80: {
     81:   PC_Factor      *dir = (PC_Factor*)pc->data;
     83:   PetscBool      flg;

     86:   if (!pc->setupcalled) {
     87:     PetscFree(dir->ordering);
     88:     PetscStrallocpy(ordering,(char**)&dir->ordering);
     89:   } else {
     90:     PetscStrcmp(dir->ordering,ordering,&flg);
     91:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change ordering after use");
     92:   }
     93:   return(0);
     94: }

     96: PetscErrorCode  PCFactorGetLevels_Factor(PC pc,PetscInt *levels)
     97: {
     98:   PC_Factor      *ilu = (PC_Factor*)pc->data;

    101:   *levels = ilu->info.levels;
    102:   return(0);
    103: }

    105: PetscErrorCode  PCFactorGetZeroPivot_Factor(PC pc,PetscReal *pivot)
    106: {
    107:   PC_Factor      *ilu = (PC_Factor*)pc->data;

    110:   *pivot = ilu->info.zeropivot;
    111:   return(0);
    112: }

    114: PetscErrorCode  PCFactorGetShiftAmount_Factor(PC pc,PetscReal *shift)
    115: {
    116:   PC_Factor      *ilu = (PC_Factor*)pc->data;

    119:   *shift = ilu->info.shiftamount;
    120:   return(0);
    121: }

    123: PetscErrorCode  PCFactorGetShiftType_Factor(PC pc,MatFactorShiftType *type)
    124: {
    125:   PC_Factor      *ilu = (PC_Factor*)pc->data;

    128:   *type = (MatFactorShiftType) (int) ilu->info.shifttype;
    129:   return(0);
    130: }

    132: PetscErrorCode  PCFactorSetLevels_Factor(PC pc,PetscInt levels)
    133: {
    134:   PC_Factor      *ilu = (PC_Factor*)pc->data;

    138:   if (!pc->setupcalled) ilu->info.levels = levels;
    139:   else if (ilu->info.levels != levels) {
    140:     (*pc->ops->reset)(pc); /* remove previous factored matrices */
    141:     pc->setupcalled  = 0; /* force a complete rebuild of preconditioner factored matrices */
    142:     ilu->info.levels = levels;
    143:   } else if (ilu->info.usedt) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change levels after use with ILUdt");
    144:   return(0);
    145: }

    147: PetscErrorCode  PCFactorSetAllowDiagonalFill_Factor(PC pc,PetscBool flg)
    148: {
    149:   PC_Factor *dir = (PC_Factor*)pc->data;

    152:   dir->info.diagonal_fill = (PetscReal) flg;
    153:   return(0);
    154: }

    156: PetscErrorCode  PCFactorGetAllowDiagonalFill_Factor(PC pc,PetscBool *flg)
    157: {
    158:   PC_Factor *dir = (PC_Factor*)pc->data;

    161:   *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
    162:   return(0);
    163: }

    165: /* ------------------------------------------------------------------------------------------*/

    167: PetscErrorCode  PCFactorSetPivotInBlocks_Factor(PC pc,PetscBool pivot)
    168: {
    169:   PC_Factor *dir = (PC_Factor*)pc->data;

    172:   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
    173:   return(0);
    174: }

    176: PetscErrorCode  PCFactorGetMatrix_Factor(PC pc,Mat *mat)
    177: {
    178:   PC_Factor *ilu = (PC_Factor*)pc->data;

    181:   if (!ilu->fact) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
    182:   *mat = ilu->fact;
    183:   return(0);
    184: }

    186: PetscErrorCode  PCFactorSetMatSolverType_Factor(PC pc,MatSolverType stype)
    187: {
    189:   PC_Factor      *lu = (PC_Factor*)pc->data;

    192:   if (lu->fact) {
    193:     MatSolverType ltype;
    194:     PetscBool     flg;
    195:     MatFactorGetSolverType(lu->fact,&ltype);
    196:     PetscStrcmp(stype,ltype,&flg);
    197:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change solver matrix package after PC has been setup or used");
    198:   }

    200:   PetscFree(lu->solvertype);
    201:   PetscStrallocpy(stype,&lu->solvertype);
    202:   return(0);
    203: }

    205: PetscErrorCode  PCFactorGetMatSolverType_Factor(PC pc,MatSolverType *stype)
    206: {
    207:   PC_Factor *lu = (PC_Factor*)pc->data;

    210:   *stype = lu->solvertype;
    211:   return(0);
    212: }

    214: PetscErrorCode  PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol)
    215: {
    216:   PC_Factor *dir = (PC_Factor*)pc->data;

    219:   if (dtcol < 0.0 || dtcol > 1.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Column pivot tolerance is %g must be between 0 and 1",(double)dtcol);
    220:   dir->info.dtcol = dtcol;
    221:   return(0);
    222: }

    224: PetscErrorCode  PCSetFromOptions_Factor(PetscOptionItems *PetscOptionsObject,PC pc)
    225: {
    226:   PC_Factor         *factor = (PC_Factor*)pc->data;
    227:   PetscErrorCode    ierr;
    228:   PetscBool         flg,set;
    229:   char              tname[256], solvertype[64];
    230:   PetscFunctionList ordlist;
    231:   PetscEnum         etmp;
    232:   PetscBool         inplace;

    235:   PCFactorGetUseInPlace(pc,&inplace);
    236:   PetscOptionsBool("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",inplace,&flg,&set);
    237:   if (set) {
    238:     PCFactorSetUseInPlace(pc,flg);
    239:   }
    240:   PetscOptionsReal("-pc_factor_fill","Expected non-zeros in factored matrix","PCFactorSetFill",((PC_Factor*)factor)->info.fill,&((PC_Factor*)factor)->info.fill,NULL);

    242:   PetscOptionsEnum("-pc_factor_shift_type","Type of shift to add to diagonal","PCFactorSetShiftType",MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);
    243:   if (flg) {
    244:     PCFactorSetShiftType(pc,(MatFactorShiftType)etmp);
    245:   }
    246:   PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,NULL);

    248:   PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,NULL);
    249:   PetscOptionsReal("-pc_factor_column_pivot","Column pivot tolerance (used only for some factorization)","PCFactorSetColumnPivot",((PC_Factor*)factor)->info.dtcol,&((PC_Factor*)factor)->info.dtcol,&flg);

    251:   PetscOptionsBool("-pc_factor_pivot_in_blocks","Pivot inside matrix dense blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",((PC_Factor*)factor)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
    252:   if (set) {
    253:     PCFactorSetPivotInBlocks(pc,flg);
    254:   }

    256:   PetscOptionsBool("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",PETSC_FALSE,&flg,&set);
    257:   if (set) {
    258:     PCFactorSetReuseFill(pc,flg);
    259:   }
    260:   PetscOptionsBool("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",PETSC_FALSE,&flg,&set);
    261:   if (set) {
    262:     PCFactorSetReuseOrdering(pc,flg);
    263:   }

    265:   MatGetOrderingList(&ordlist);
    266:   PetscOptionsFList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in factored matrix","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)factor)->ordering,tname,sizeof(tname),&flg);
    267:   if (flg) {
    268:     PCFactorSetMatOrderingType(pc,tname);
    269:   }

    271:   /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
    272:   PetscOptionsDeprecated("-pc_factor_mat_solver_package","-pc_factor_mat_solver_type","3.9",NULL);
    273:   PetscOptionsString("-pc_factor_mat_solver_type","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,sizeof(solvertype),&flg);
    274:   if (flg) {
    275:     PCFactorSetMatSolverType(pc,solvertype);
    276:   }
    277:   return(0);
    278: }

    280: PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer)
    281: {
    282:   PC_Factor       *factor = (PC_Factor*)pc->data;
    283:   PetscErrorCode  ierr;
    284:   PetscBool       isstring,iascii,flg;
    285:   MatOrderingType ordering;

    288:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
    289:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    290:   if (iascii) {
    291:     if (factor->inplace) {
    292:       PetscViewerASCIIPrintf(viewer,"  in-place factorization\n");
    293:     } else {
    294:       PetscViewerASCIIPrintf(viewer,"  out-of-place factorization\n");
    295:     }

    297:     if (factor->reusefill)     {PetscViewerASCIIPrintf(viewer,"  Reusing fill from past factorization\n");}
    298:     if (factor->reuseordering) {PetscViewerASCIIPrintf(viewer,"  Reusing reordering from past factorization\n");}
    299:     if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
    300:       if (factor->info.dt > 0) {
    301:         PetscViewerASCIIPrintf(viewer,"  drop tolerance %g\n",(double)factor->info.dt);
    302:         PetscViewerASCIIPrintf(viewer,"  max nonzeros per row %D\n",factor->info.dtcount);
    303:         PetscViewerASCIIPrintf(viewer,"  column permutation tolerance %g\n",(double)factor->info.dtcol);
    304:       } else if (factor->info.levels == 1) {
    305:         PetscViewerASCIIPrintf(viewer,"  %D level of fill\n",(PetscInt)factor->info.levels);
    306:       } else {
    307:         PetscViewerASCIIPrintf(viewer,"  %D levels of fill\n",(PetscInt)factor->info.levels);
    308:       }
    309:     }

    311:     PetscViewerASCIIPrintf(viewer,"  tolerance for zero pivot %g\n",(double)factor->info.zeropivot);
    312:     if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */
    313:       PetscViewerASCIIPrintf(viewer,"  using %s [%s]\n",MatFactorShiftTypesDetail[(int)factor->info.shifttype],MatFactorShiftTypes[(int)factor->info.shifttype]);
    314:     }

    316:     PetscStrcmp(factor->ordering,MATORDERINGNATURAL_OR_ND,&flg);
    317:     if (flg) {
    318:       PetscBool isseqsbaij;
    319:       PetscObjectTypeCompareAny((PetscObject)pc->pmat,&isseqsbaij,MATSEQSBAIJ,MATSEQBAIJ,NULL);
    320:       if (isseqsbaij) {
    321:         ordering = MATORDERINGNATURAL;
    322:       } else {
    323:         ordering = MATORDERINGND;
    324:       }
    325:     } else {
    326:       ordering = factor->ordering;
    327:     }
    328:     if (!factor->fact) {
    329:       PetscViewerASCIIPrintf(viewer,"  matrix ordering: %s (may be overridden during setup)\n",ordering);
    330:     } else {
    331:       PetscBool useordering;
    332:       MatInfo   info;
    333:       MatFactorGetUseOrdering(factor->fact,&useordering);
    334:       if (!useordering) ordering = "external";
    335:       PetscViewerASCIIPrintf(viewer,"  matrix ordering: %s\n",ordering);
    336:       MatGetInfo(factor->fact,MAT_LOCAL,&info);
    337:       PetscViewerASCIIPrintf(viewer,"  factor fill ratio given %g, needed %g\n",(double)info.fill_ratio_given,(double)info.fill_ratio_needed);
    338:       PetscViewerASCIIPrintf(viewer,"    Factored matrix follows:\n");
    339:       PetscViewerASCIIPushTab(viewer);
    340:       PetscViewerASCIIPushTab(viewer);
    341:       PetscViewerASCIIPushTab(viewer);
    342:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
    343:       MatView(factor->fact,viewer);
    344:       PetscViewerPopFormat(viewer);
    345:       PetscViewerASCIIPopTab(viewer);
    346:       PetscViewerASCIIPopTab(viewer);
    347:       PetscViewerASCIIPopTab(viewer);
    348:     }

    350:   } else if (isstring) {
    351:     MatFactorType t;
    352:     MatGetFactorType(factor->fact,&t);
    353:     if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) {
    354:       PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);
    355:     }
    356:   }
    357:   return(0);
    358: }