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