Actual source code: factimpl.c

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

      3: PetscErrorCode PCFactorSetUpMatSolverType_Factor(PC pc)
      4: {
      5:   PC_Factor *icc = (PC_Factor *)pc->data;

      7:   PetscFunctionBegin;
      8:   PetscCheck(pc->pmat, 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()");
      9:   if (!((PC_Factor *)icc)->fact) PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, ((PC_Factor *)icc)->factortype, &((PC_Factor *)icc)->fact));
     10:   PetscCheck(((PC_Factor *)icc)->fact, PetscObjectComm((PetscObject)pc->pmat), PETSC_ERR_SUP, "MatFactor type %s not supported by this matrix instance of type %s and solver type %s. Rerun with \"-info :mat | grep MatGetFactor_\" for additional information.",
     11:              MatFactorTypes[((PC_Factor *)icc)->factortype], ((PetscObject)pc->pmat)->type_name, ((PC_Factor *)icc)->solvertype);
     12:   PetscFunctionReturn(PETSC_SUCCESS);
     13: }

     15: PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc, PetscReal z)
     16: {
     17:   PC_Factor *ilu = (PC_Factor *)pc->data;

     19:   PetscFunctionBegin;
     20:   ilu->info.zeropivot = z;
     21:   PetscFunctionReturn(PETSC_SUCCESS);
     22: }

     24: PetscErrorCode PCFactorSetShiftType_Factor(PC pc, MatFactorShiftType shifttype)
     25: {
     26:   PC_Factor *dir = (PC_Factor *)pc->data;

     28:   PetscFunctionBegin;
     29:   if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
     30:   else {
     31:     dir->info.shifttype = (PetscReal)shifttype;
     32:     if ((shifttype == MAT_SHIFT_NONZERO || shifttype == MAT_SHIFT_INBLOCKS) && dir->info.shiftamount == 0.0) { dir->info.shiftamount = 100.0 * PETSC_MACHINE_EPSILON; /* set default amount if user has not called PCFactorSetShiftAmount() yet */ }
     33:   }
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc, PetscReal shiftamount)
     38: {
     39:   PC_Factor *dir = (PC_Factor *)pc->data;

     41:   PetscFunctionBegin;
     42:   if (shiftamount == (PetscReal)PETSC_DECIDE) dir->info.shiftamount = 100.0 * PETSC_MACHINE_EPSILON;
     43:   else dir->info.shiftamount = shiftamount;
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
     48: {
     49:   PC_Factor *ilu = (PC_Factor *)pc->data;

     51:   PetscFunctionBegin;
     52:   if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor *)ilu)->info.dt != dt || ((PC_Factor *)ilu)->info.dtcol != dtcol || ((PC_Factor *)ilu)->info.dtcount != dtcount)) {
     53:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change tolerance after use");
     54:   }
     55:   ilu->info.usedt   = PETSC_TRUE;
     56:   ilu->info.dt      = dt;
     57:   ilu->info.dtcol   = dtcol;
     58:   ilu->info.dtcount = dtcount;
     59:   ilu->info.fill    = PETSC_DEFAULT;
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: PetscErrorCode PCFactorSetFill_Factor(PC pc, PetscReal fill)
     64: {
     65:   PC_Factor *dir = (PC_Factor *)pc->data;

     67:   PetscFunctionBegin;
     68:   dir->info.fill = fill;
     69:   PetscFunctionReturn(PETSC_SUCCESS);
     70: }

     72: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc, MatOrderingType ordering)
     73: {
     74:   PC_Factor *dir = (PC_Factor *)pc->data;
     75:   PetscBool  flg;

     77:   PetscFunctionBegin;
     78:   if (!pc->setupcalled) {
     79:     PetscCall(PetscFree(dir->ordering));
     80:     PetscCall(PetscStrallocpy(ordering, (char **)&dir->ordering));
     81:   } else {
     82:     PetscCall(PetscStrcmp(dir->ordering, ordering, &flg));
     83:     PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change ordering after use");
     84:   }
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: PetscErrorCode PCFactorGetLevels_Factor(PC pc, PetscInt *levels)
     89: {
     90:   PC_Factor *ilu = (PC_Factor *)pc->data;

     92:   PetscFunctionBegin;
     93:   *levels = ilu->info.levels;
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: PetscErrorCode PCFactorGetZeroPivot_Factor(PC pc, PetscReal *pivot)
     98: {
     99:   PC_Factor *ilu = (PC_Factor *)pc->data;

    101:   PetscFunctionBegin;
    102:   *pivot = ilu->info.zeropivot;
    103:   PetscFunctionReturn(PETSC_SUCCESS);
    104: }

    106: PetscErrorCode PCFactorGetShiftAmount_Factor(PC pc, PetscReal *shift)
    107: {
    108:   PC_Factor *ilu = (PC_Factor *)pc->data;

    110:   PetscFunctionBegin;
    111:   *shift = ilu->info.shiftamount;
    112:   PetscFunctionReturn(PETSC_SUCCESS);
    113: }

    115: PetscErrorCode PCFactorGetShiftType_Factor(PC pc, MatFactorShiftType *type)
    116: {
    117:   PC_Factor *ilu = (PC_Factor *)pc->data;

    119:   PetscFunctionBegin;
    120:   *type = (MatFactorShiftType)(int)ilu->info.shifttype;
    121:   PetscFunctionReturn(PETSC_SUCCESS);
    122: }

    124: PetscErrorCode PCFactorSetLevels_Factor(PC pc, PetscInt levels)
    125: {
    126:   PC_Factor *ilu = (PC_Factor *)pc->data;

    128:   PetscFunctionBegin;
    129:   if (!pc->setupcalled) ilu->info.levels = levels;
    130:   else if (ilu->info.levels != levels) {
    131:     PetscUseTypeMethod(pc, reset); /* remove previous factored matrices */
    132:     pc->setupcalled  = 0;          /* force a complete rebuild of preconditioner factored matrices */
    133:     ilu->info.levels = levels;
    134:   } else PetscCheck(!ilu->info.usedt, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change levels after use with ILUdt");
    135:   PetscFunctionReturn(PETSC_SUCCESS);
    136: }

    138: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc, PetscBool flg)
    139: {
    140:   PC_Factor *dir = (PC_Factor *)pc->data;

    142:   PetscFunctionBegin;
    143:   dir->info.diagonal_fill = (PetscReal)flg;
    144:   PetscFunctionReturn(PETSC_SUCCESS);
    145: }

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

    151:   PetscFunctionBegin;
    152:   *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
    153:   PetscFunctionReturn(PETSC_SUCCESS);
    154: }

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

    160:   PetscFunctionBegin;
    161:   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
    162:   PetscFunctionReturn(PETSC_SUCCESS);
    163: }

    165: PetscErrorCode PCFactorGetMatrix_Factor(PC pc, Mat *mat)
    166: {
    167:   PC_Factor *ilu = (PC_Factor *)pc->data;

    169:   PetscFunctionBegin;
    170:   PetscCheck(ilu->fact, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
    171:   *mat = ilu->fact;
    172:   PetscFunctionReturn(PETSC_SUCCESS);
    173: }

    175: /* allow access to preallocation information */
    176: #include <petsc/private/matimpl.h>

    178: PetscErrorCode PCFactorSetMatSolverType_Factor(PC pc, MatSolverType stype)
    179: {
    180:   PC_Factor *lu = (PC_Factor *)pc->data;

    182:   PetscFunctionBegin;
    183:   if (lu->fact && lu->fact->assembled) {
    184:     MatSolverType ltype;
    185:     PetscBool     flg;
    186:     PetscCall(MatFactorGetSolverType(lu->fact, &ltype));
    187:     PetscCall(PetscStrcmp(stype, ltype, &flg));
    188:     PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change solver matrix package from %s to %s after PC has been setup or used", ltype, stype);
    189:   }

    191:   PetscCall(PetscFree(lu->solvertype));
    192:   PetscCall(PetscStrallocpy(stype, &lu->solvertype));
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }

    196: PetscErrorCode PCFactorGetMatSolverType_Factor(PC pc, MatSolverType *stype)
    197: {
    198:   PC_Factor *lu = (PC_Factor *)pc->data;

    200:   PetscFunctionBegin;
    201:   *stype = lu->solvertype;
    202:   PetscFunctionReturn(PETSC_SUCCESS);
    203: }

    205: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc, PetscReal dtcol)
    206: {
    207:   PC_Factor *dir = (PC_Factor *)pc->data;

    209:   PetscFunctionBegin;
    210:   PetscCheck(dtcol >= 0.0 && dtcol <= 1.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Column pivot tolerance is %g must be between 0 and 1", (double)dtcol);
    211:   dir->info.dtcol = dtcol;
    212:   PetscFunctionReturn(PETSC_SUCCESS);
    213: }

    215: PetscErrorCode PCSetFromOptions_Factor(PC pc, PetscOptionItems *PetscOptionsObject)
    216: {
    217:   PC_Factor        *factor = (PC_Factor *)pc->data;
    218:   PetscBool         flg, set;
    219:   char              tname[256], solvertype[64];
    220:   PetscFunctionList ordlist;
    221:   PetscEnum         etmp;
    222:   PetscBool         inplace;

    224:   PetscFunctionBegin;
    225:   PetscCall(PCFactorGetUseInPlace(pc, &inplace));
    226:   PetscCall(PetscOptionsBool("-pc_factor_in_place", "Form factored matrix in the same memory as the matrix", "PCFactorSetUseInPlace", inplace, &flg, &set));
    227:   if (set) PetscCall(PCFactorSetUseInPlace(pc, flg));
    228:   PetscCall(PetscOptionsReal("-pc_factor_fill", "Expected non-zeros in factored matrix", "PCFactorSetFill", ((PC_Factor *)factor)->info.fill, &((PC_Factor *)factor)->info.fill, NULL));

    230:   PetscCall(PetscOptionsEnum("-pc_factor_shift_type", "Type of shift to add to diagonal", "PCFactorSetShiftType", MatFactorShiftTypes, (PetscEnum)(int)((PC_Factor *)factor)->info.shifttype, &etmp, &flg));
    231:   if (flg) PetscCall(PCFactorSetShiftType(pc, (MatFactorShiftType)etmp));
    232:   PetscCall(PetscOptionsReal("-pc_factor_shift_amount", "Shift added to diagonal", "PCFactorSetShiftAmount", ((PC_Factor *)factor)->info.shiftamount, &((PC_Factor *)factor)->info.shiftamount, NULL));

    234:   PetscCall(PetscOptionsReal("-pc_factor_zeropivot", "Pivot is considered zero if less than", "PCFactorSetZeroPivot", ((PC_Factor *)factor)->info.zeropivot, &((PC_Factor *)factor)->info.zeropivot, NULL));
    235:   PetscCall(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));

    237:   PetscCall(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));
    238:   if (set) PetscCall(PCFactorSetPivotInBlocks(pc, flg));

    240:   PetscCall(PetscOptionsBool("-pc_factor_reuse_fill", "Use fill from previous factorization", "PCFactorSetReuseFill", PETSC_FALSE, &flg, &set));
    241:   if (set) PetscCall(PCFactorSetReuseFill(pc, flg));
    242:   PetscCall(PetscOptionsBool("-pc_factor_reuse_ordering", "Reuse ordering from previous factorization", "PCFactorSetReuseOrdering", PETSC_FALSE, &flg, &set));
    243:   if (set) PetscCall(PCFactorSetReuseOrdering(pc, flg));

    245:   PetscCall(PetscOptionsDeprecated("-pc_factor_mat_solver_package", "-pc_factor_mat_solver_type", "3.9", NULL));
    246:   PetscCall(PetscOptionsString("-pc_factor_mat_solver_type", "Specific direct solver to use", "MatGetFactor", ((PC_Factor *)factor)->solvertype, solvertype, sizeof(solvertype), &flg));
    247:   if (flg) PetscCall(PCFactorSetMatSolverType(pc, solvertype));
    248:   PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
    249:   PetscCall(MatGetOrderingList(&ordlist));
    250:   PetscCall(PetscOptionsFList("-pc_factor_mat_ordering_type", "Reordering to reduce nonzeros in factored matrix", "PCFactorSetMatOrderingType", ordlist, ((PC_Factor *)factor)->ordering, tname, sizeof(tname), &flg));
    251:   if (flg) PetscCall(PCFactorSetMatOrderingType(pc, tname));
    252:   PetscFunctionReturn(PETSC_SUCCESS);
    253: }

    255: PetscErrorCode PCView_Factor(PC pc, PetscViewer viewer)
    256: {
    257:   PC_Factor        *factor = (PC_Factor *)pc->data;
    258:   PetscBool         isstring, iascii, canuseordering;
    259:   MatInfo           info;
    260:   MatOrderingType   ordering;
    261:   PetscViewerFormat format;

    263:   PetscFunctionBegin;
    264:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
    265:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    266:   if (iascii) {
    267:     if (factor->inplace) {
    268:       PetscCall(PetscViewerASCIIPrintf(viewer, "  in-place factorization\n"));
    269:     } else {
    270:       PetscCall(PetscViewerASCIIPrintf(viewer, "  out-of-place factorization\n"));
    271:     }

    273:     if (factor->reusefill) PetscCall(PetscViewerASCIIPrintf(viewer, "  Reusing fill from past factorization\n"));
    274:     if (factor->reuseordering) PetscCall(PetscViewerASCIIPrintf(viewer, "  Reusing reordering from past factorization\n"));
    275:     if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
    276:       if (factor->info.dt > 0) {
    277:         PetscCall(PetscViewerASCIIPrintf(viewer, "  drop tolerance %g\n", (double)factor->info.dt));
    278:         PetscCall(PetscViewerASCIIPrintf(viewer, "  max nonzeros per row %" PetscInt_FMT "\n", (PetscInt)factor->info.dtcount));
    279:         PetscCall(PetscViewerASCIIPrintf(viewer, "  column permutation tolerance %g\n", (double)factor->info.dtcol));
    280:       } else if (factor->info.levels == 1) {
    281:         PetscCall(PetscViewerASCIIPrintf(viewer, "  %" PetscInt_FMT " level of fill\n", (PetscInt)factor->info.levels));
    282:       } else {
    283:         PetscCall(PetscViewerASCIIPrintf(viewer, "  %" PetscInt_FMT " levels of fill\n", (PetscInt)factor->info.levels));
    284:       }
    285:     }

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

    292:     if (factor->fact) {
    293:       PetscCall(MatFactorGetCanUseOrdering(factor->fact, &canuseordering));
    294:       if (!canuseordering) ordering = MATORDERINGEXTERNAL;
    295:       else ordering = factor->ordering;
    296:       PetscCall(PetscViewerASCIIPrintf(viewer, "  matrix ordering: %s\n", ordering));
    297:       if (!factor->fact->assembled) {
    298:         PetscCall(PetscViewerASCIIPrintf(viewer, "  matrix solver type: %s\n", factor->fact->solvertype));
    299:         PetscCall(PetscViewerASCIIPrintf(viewer, "  matrix not yet factored; no additional information available\n"));
    300:       } else {
    301:         PetscCall(MatGetInfo(factor->fact, MAT_LOCAL, &info));
    302:         PetscCall(PetscViewerASCIIPrintf(viewer, "  factor fill ratio given %g, needed %g\n", (double)info.fill_ratio_given, (double)info.fill_ratio_needed));
    303:         PetscCall(PetscViewerASCIIPrintf(viewer, "    Factored matrix follows:\n"));
    304:         PetscCall(PetscViewerASCIIPushTab(viewer));
    305:         PetscCall(PetscViewerASCIIPushTab(viewer));
    306:         PetscCall(PetscViewerASCIIPushTab(viewer));
    307:         PetscCall(PetscViewerGetFormat(viewer, &format));
    308:         PetscCall(PetscViewerPushFormat(viewer, format != PETSC_VIEWER_ASCII_INFO_DETAIL ? PETSC_VIEWER_ASCII_INFO : PETSC_VIEWER_ASCII_INFO_DETAIL));
    309:         PetscCall(MatView(factor->fact, viewer));
    310:         PetscCall(PetscViewerPopFormat(viewer));
    311:         PetscCall(PetscViewerASCIIPopTab(viewer));
    312:         PetscCall(PetscViewerASCIIPopTab(viewer));
    313:         PetscCall(PetscViewerASCIIPopTab(viewer));
    314:       }
    315:     }

    317:   } else if (isstring) {
    318:     MatFactorType t;
    319:     PetscCall(MatGetFactorType(factor->fact, &t));
    320:     if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) PetscCall(PetscViewerStringSPrintf(viewer, " lvls=%" PetscInt_FMT ",order=%s", (PetscInt)factor->info.levels, factor->ordering));
    321:   }
    322:   PetscFunctionReturn(PETSC_SUCCESS);
    323: }