1 2 #include <../src/ksp/pc/impls/factor/factor.h> /*I "petscpc.h" I*/ 3 4 PetscErrorCode PCFactorSetUpMatSolverType_Factor(PC pc) { 5 PC_Factor *icc = (PC_Factor *)pc->data; 6 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->setupcalled && !((PC_Factor *)icc)->fact) PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, ((PC_Factor *)icc)->factortype, &((PC_Factor *)icc)->fact)); 10 PetscFunctionReturn(0); 11 } 12 13 PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc, PetscReal z) { 14 PC_Factor *ilu = (PC_Factor *)pc->data; 15 16 PetscFunctionBegin; 17 ilu->info.zeropivot = z; 18 PetscFunctionReturn(0); 19 } 20 21 PetscErrorCode PCFactorSetShiftType_Factor(PC pc, MatFactorShiftType shifttype) { 22 PC_Factor *dir = (PC_Factor *)pc->data; 23 24 PetscFunctionBegin; 25 if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal)MAT_SHIFT_NONE; 26 else { 27 dir->info.shifttype = (PetscReal)shifttype; 28 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 */ } 29 } 30 PetscFunctionReturn(0); 31 } 32 33 PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc, PetscReal shiftamount) { 34 PC_Factor *dir = (PC_Factor *)pc->data; 35 36 PetscFunctionBegin; 37 if (shiftamount == (PetscReal)PETSC_DECIDE) dir->info.shiftamount = 100.0 * PETSC_MACHINE_EPSILON; 38 else dir->info.shiftamount = shiftamount; 39 PetscFunctionReturn(0); 40 } 41 42 PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount) { 43 PC_Factor *ilu = (PC_Factor *)pc->data; 44 45 PetscFunctionBegin; 46 if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor *)ilu)->info.dt != dt || ((PC_Factor *)ilu)->info.dtcol != dtcol || ((PC_Factor *)ilu)->info.dtcount != dtcount)) { 47 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change tolerance after use"); 48 } 49 ilu->info.usedt = PETSC_TRUE; 50 ilu->info.dt = dt; 51 ilu->info.dtcol = dtcol; 52 ilu->info.dtcount = dtcount; 53 ilu->info.fill = PETSC_DEFAULT; 54 PetscFunctionReturn(0); 55 } 56 57 PetscErrorCode PCFactorSetFill_Factor(PC pc, PetscReal fill) { 58 PC_Factor *dir = (PC_Factor *)pc->data; 59 60 PetscFunctionBegin; 61 dir->info.fill = fill; 62 PetscFunctionReturn(0); 63 } 64 65 PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc, MatOrderingType ordering) { 66 PC_Factor *dir = (PC_Factor *)pc->data; 67 PetscBool flg; 68 69 PetscFunctionBegin; 70 if (!pc->setupcalled) { 71 PetscCall(PetscFree(dir->ordering)); 72 PetscCall(PetscStrallocpy(ordering, (char **)&dir->ordering)); 73 } else { 74 PetscCall(PetscStrcmp(dir->ordering, ordering, &flg)); 75 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change ordering after use"); 76 } 77 PetscFunctionReturn(0); 78 } 79 80 PetscErrorCode PCFactorGetLevels_Factor(PC pc, PetscInt *levels) { 81 PC_Factor *ilu = (PC_Factor *)pc->data; 82 83 PetscFunctionBegin; 84 *levels = ilu->info.levels; 85 PetscFunctionReturn(0); 86 } 87 88 PetscErrorCode PCFactorGetZeroPivot_Factor(PC pc, PetscReal *pivot) { 89 PC_Factor *ilu = (PC_Factor *)pc->data; 90 91 PetscFunctionBegin; 92 *pivot = ilu->info.zeropivot; 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode PCFactorGetShiftAmount_Factor(PC pc, PetscReal *shift) { 97 PC_Factor *ilu = (PC_Factor *)pc->data; 98 99 PetscFunctionBegin; 100 *shift = ilu->info.shiftamount; 101 PetscFunctionReturn(0); 102 } 103 104 PetscErrorCode PCFactorGetShiftType_Factor(PC pc, MatFactorShiftType *type) { 105 PC_Factor *ilu = (PC_Factor *)pc->data; 106 107 PetscFunctionBegin; 108 *type = (MatFactorShiftType)(int)ilu->info.shifttype; 109 PetscFunctionReturn(0); 110 } 111 112 PetscErrorCode PCFactorSetLevels_Factor(PC pc, PetscInt levels) { 113 PC_Factor *ilu = (PC_Factor *)pc->data; 114 115 PetscFunctionBegin; 116 if (!pc->setupcalled) ilu->info.levels = levels; 117 else if (ilu->info.levels != levels) { 118 PetscUseTypeMethod(pc, reset); /* remove previous factored matrices */ 119 pc->setupcalled = 0; /* force a complete rebuild of preconditioner factored matrices */ 120 ilu->info.levels = levels; 121 } else PetscCheck(!ilu->info.usedt, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change levels after use with ILUdt"); 122 PetscFunctionReturn(0); 123 } 124 125 PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc, PetscBool flg) { 126 PC_Factor *dir = (PC_Factor *)pc->data; 127 128 PetscFunctionBegin; 129 dir->info.diagonal_fill = (PetscReal)flg; 130 PetscFunctionReturn(0); 131 } 132 133 PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc, PetscBool *flg) { 134 PC_Factor *dir = (PC_Factor *)pc->data; 135 136 PetscFunctionBegin; 137 *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE; 138 PetscFunctionReturn(0); 139 } 140 141 PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc, PetscBool pivot) { 142 PC_Factor *dir = (PC_Factor *)pc->data; 143 144 PetscFunctionBegin; 145 dir->info.pivotinblocks = pivot ? 1.0 : 0.0; 146 PetscFunctionReturn(0); 147 } 148 149 PetscErrorCode PCFactorGetMatrix_Factor(PC pc, Mat *mat) { 150 PC_Factor *ilu = (PC_Factor *)pc->data; 151 152 PetscFunctionBegin; 153 PetscCheck(ilu->fact, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Matrix not yet factored; call after KSPSetUp() or PCSetUp()"); 154 *mat = ilu->fact; 155 PetscFunctionReturn(0); 156 } 157 158 /* allow access to preallocation information */ 159 #include <petsc/private/matimpl.h> 160 161 PetscErrorCode PCFactorSetMatSolverType_Factor(PC pc, MatSolverType stype) { 162 PC_Factor *lu = (PC_Factor *)pc->data; 163 164 PetscFunctionBegin; 165 if (lu->fact && lu->fact->assembled) { 166 MatSolverType ltype; 167 PetscBool flg; 168 PetscCall(MatFactorGetSolverType(lu->fact, <ype)); 169 PetscCall(PetscStrcmp(stype, ltype, &flg)); 170 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); 171 } 172 173 PetscCall(PetscFree(lu->solvertype)); 174 PetscCall(PetscStrallocpy(stype, &lu->solvertype)); 175 PetscFunctionReturn(0); 176 } 177 178 PetscErrorCode PCFactorGetMatSolverType_Factor(PC pc, MatSolverType *stype) { 179 PC_Factor *lu = (PC_Factor *)pc->data; 180 181 PetscFunctionBegin; 182 *stype = lu->solvertype; 183 PetscFunctionReturn(0); 184 } 185 186 PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc, PetscReal dtcol) { 187 PC_Factor *dir = (PC_Factor *)pc->data; 188 189 PetscFunctionBegin; 190 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); 191 dir->info.dtcol = dtcol; 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode PCSetFromOptions_Factor(PC pc, PetscOptionItems *PetscOptionsObject) { 196 PC_Factor *factor = (PC_Factor *)pc->data; 197 PetscBool flg, set; 198 char tname[256], solvertype[64]; 199 PetscFunctionList ordlist; 200 PetscEnum etmp; 201 PetscBool inplace; 202 203 PetscFunctionBegin; 204 PetscCall(PCFactorGetUseInPlace(pc, &inplace)); 205 PetscCall(PetscOptionsBool("-pc_factor_in_place", "Form factored matrix in the same memory as the matrix", "PCFactorSetUseInPlace", inplace, &flg, &set)); 206 if (set) PetscCall(PCFactorSetUseInPlace(pc, flg)); 207 PetscCall(PetscOptionsReal("-pc_factor_fill", "Expected non-zeros in factored matrix", "PCFactorSetFill", ((PC_Factor *)factor)->info.fill, &((PC_Factor *)factor)->info.fill, NULL)); 208 209 PetscCall(PetscOptionsEnum("-pc_factor_shift_type", "Type of shift to add to diagonal", "PCFactorSetShiftType", MatFactorShiftTypes, (PetscEnum)(int)((PC_Factor *)factor)->info.shifttype, &etmp, &flg)); 210 if (flg) PetscCall(PCFactorSetShiftType(pc, (MatFactorShiftType)etmp)); 211 PetscCall(PetscOptionsReal("-pc_factor_shift_amount", "Shift added to diagonal", "PCFactorSetShiftAmount", ((PC_Factor *)factor)->info.shiftamount, &((PC_Factor *)factor)->info.shiftamount, NULL)); 212 213 PetscCall(PetscOptionsReal("-pc_factor_zeropivot", "Pivot is considered zero if less than", "PCFactorSetZeroPivot", ((PC_Factor *)factor)->info.zeropivot, &((PC_Factor *)factor)->info.zeropivot, NULL)); 214 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)); 215 216 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)); 217 if (set) PetscCall(PCFactorSetPivotInBlocks(pc, flg)); 218 219 PetscCall(PetscOptionsBool("-pc_factor_reuse_fill", "Use fill from previous factorization", "PCFactorSetReuseFill", PETSC_FALSE, &flg, &set)); 220 if (set) PetscCall(PCFactorSetReuseFill(pc, flg)); 221 PetscCall(PetscOptionsBool("-pc_factor_reuse_ordering", "Reuse ordering from previous factorization", "PCFactorSetReuseOrdering", PETSC_FALSE, &flg, &set)); 222 if (set) PetscCall(PCFactorSetReuseOrdering(pc, flg)); 223 224 PetscCall(PetscOptionsDeprecated("-pc_factor_mat_solver_package", "-pc_factor_mat_solver_type", "3.9", NULL)); 225 PetscCall(PetscOptionsString("-pc_factor_mat_solver_type", "Specific direct solver to use", "MatGetFactor", ((PC_Factor *)factor)->solvertype, solvertype, sizeof(solvertype), &flg)); 226 if (flg) PetscCall(PCFactorSetMatSolverType(pc, solvertype)); 227 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 228 PetscCall(MatGetOrderingList(&ordlist)); 229 PetscCall(PetscOptionsFList("-pc_factor_mat_ordering_type", "Reordering to reduce nonzeros in factored matrix", "PCFactorSetMatOrderingType", ordlist, ((PC_Factor *)factor)->ordering, tname, sizeof(tname), &flg)); 230 if (flg) PetscCall(PCFactorSetMatOrderingType(pc, tname)); 231 PetscFunctionReturn(0); 232 } 233 234 PetscErrorCode PCView_Factor(PC pc, PetscViewer viewer) { 235 PC_Factor *factor = (PC_Factor *)pc->data; 236 PetscBool isstring, iascii, canuseordering; 237 MatInfo info; 238 MatOrderingType ordering; 239 PetscViewerFormat format; 240 241 PetscFunctionBegin; 242 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring)); 243 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 244 if (iascii) { 245 if (factor->inplace) { 246 PetscCall(PetscViewerASCIIPrintf(viewer, " in-place factorization\n")); 247 } else { 248 PetscCall(PetscViewerASCIIPrintf(viewer, " out-of-place factorization\n")); 249 } 250 251 if (factor->reusefill) PetscCall(PetscViewerASCIIPrintf(viewer, " Reusing fill from past factorization\n")); 252 if (factor->reuseordering) PetscCall(PetscViewerASCIIPrintf(viewer, " Reusing reordering from past factorization\n")); 253 if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) { 254 if (factor->info.dt > 0) { 255 PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)factor->info.dt)); 256 PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzeros per row %" PetscInt_FMT "\n", (PetscInt)factor->info.dtcount)); 257 PetscCall(PetscViewerASCIIPrintf(viewer, " column permutation tolerance %g\n", (double)factor->info.dtcol)); 258 } else if (factor->info.levels == 1) { 259 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " level of fill\n", (PetscInt)factor->info.levels)); 260 } else { 261 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " levels of fill\n", (PetscInt)factor->info.levels)); 262 } 263 } 264 265 PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance for zero pivot %g\n", (double)factor->info.zeropivot)); 266 if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */ 267 PetscCall(PetscViewerASCIIPrintf(viewer, " using %s [%s]\n", MatFactorShiftTypesDetail[(int)factor->info.shifttype], MatFactorShiftTypes[(int)factor->info.shifttype])); 268 } 269 270 if (factor->fact) { 271 PetscCall(MatFactorGetCanUseOrdering(factor->fact, &canuseordering)); 272 if (!canuseordering) ordering = MATORDERINGEXTERNAL; 273 else ordering = factor->ordering; 274 PetscCall(PetscViewerASCIIPrintf(viewer, " matrix ordering: %s\n", ordering)); 275 if (!factor->fact->assembled) { 276 PetscCall(PetscViewerASCIIPrintf(viewer, " matrix solver type: %s\n", factor->fact->solvertype)); 277 PetscCall(PetscViewerASCIIPrintf(viewer, " matrix not yet factored; no additional information available\n")); 278 } else { 279 PetscCall(MatGetInfo(factor->fact, MAT_LOCAL, &info)); 280 PetscCall(PetscViewerASCIIPrintf(viewer, " factor fill ratio given %g, needed %g\n", (double)info.fill_ratio_given, (double)info.fill_ratio_needed)); 281 PetscCall(PetscViewerASCIIPrintf(viewer, " Factored matrix follows:\n")); 282 PetscCall(PetscViewerASCIIPushTab(viewer)); 283 PetscCall(PetscViewerASCIIPushTab(viewer)); 284 PetscCall(PetscViewerASCIIPushTab(viewer)); 285 PetscCall(PetscViewerGetFormat(viewer, &format)); 286 PetscCall(PetscViewerPushFormat(viewer, format != PETSC_VIEWER_ASCII_INFO_DETAIL ? PETSC_VIEWER_ASCII_INFO : PETSC_VIEWER_ASCII_INFO_DETAIL)); 287 PetscCall(MatView(factor->fact, viewer)); 288 PetscCall(PetscViewerPopFormat(viewer)); 289 PetscCall(PetscViewerASCIIPopTab(viewer)); 290 PetscCall(PetscViewerASCIIPopTab(viewer)); 291 PetscCall(PetscViewerASCIIPopTab(viewer)); 292 } 293 } 294 295 } else if (isstring) { 296 MatFactorType t; 297 PetscCall(MatGetFactorType(factor->fact, &t)); 298 if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) PetscCall(PetscViewerStringSPrintf(viewer, " lvls=%" PetscInt_FMT ",order=%s", (PetscInt)factor->info.levels, factor->ordering)); 299 } 300 PetscFunctionReturn(0); 301 } 302