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