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