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