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