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