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