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