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