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