1 #define PETSCKSP_DLL 2 3 /* 4 Defines a ILU factorization preconditioner for any Mat implementation 5 */ 6 #include "../src/ksp/pc/impls/factor/ilu/ilu.h" /*I "petscpc.h" I*/ 7 8 /* ------------------------------------------------------------------------------------------*/ 9 EXTERN_C_BEGIN 10 #undef __FUNCT__ 11 #define __FUNCT__ "PCFactorSetReuseFill_ILU" 12 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseFill_ILU(PC pc,PetscTruth flag) 13 { 14 PC_ILU *lu = (PC_ILU*)pc->data; 15 16 PetscFunctionBegin; 17 lu->reusefill = flag; 18 PetscFunctionReturn(0); 19 } 20 EXTERN_C_END 21 22 EXTERN_C_BEGIN 23 #undef __FUNCT__ 24 #define __FUNCT__ "PCFactorReorderForNonzeroDiagonal_ILU" 25 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z) 26 { 27 PC_ILU *ilu = (PC_ILU*)pc->data; 28 29 PetscFunctionBegin; 30 ilu->nonzerosalongdiagonal = PETSC_TRUE; 31 if (z == PETSC_DECIDE) { 32 ilu->nonzerosalongdiagonaltol = 1.e-10; 33 } else { 34 ilu->nonzerosalongdiagonaltol = z; 35 } 36 PetscFunctionReturn(0); 37 } 38 EXTERN_C_END 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "PCDestroy_ILU_Internal" 42 PetscErrorCode PCDestroy_ILU_Internal(PC pc) 43 { 44 PC_ILU *ilu = (PC_ILU*)pc->data; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 if (!ilu->inplace && ((PC_Factor*)ilu)->fact) {ierr = MatDestroy(((PC_Factor*)ilu)->fact);CHKERRQ(ierr);} 49 if (ilu->row && ilu->col && ilu->row != ilu->col) {ierr = ISDestroy(ilu->row);CHKERRQ(ierr);} 50 if (ilu->col) {ierr = ISDestroy(ilu->col);CHKERRQ(ierr);} 51 PetscFunctionReturn(0); 52 } 53 54 EXTERN_C_BEGIN 55 #undef __FUNCT__ 56 #define __FUNCT__ "PCFactorSetUseDropTolerance_ILU" 57 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount) 58 { 59 PC_ILU *ilu = (PC_ILU*)pc->data; 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 if (pc->setupcalled && (!ilu->usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) { 64 pc->setupcalled = 0; 65 ierr = PCDestroy_ILU_Internal(pc);CHKERRQ(ierr); 66 } 67 ilu->usedt = PETSC_TRUE; 68 ((PC_Factor*)ilu)->info.dt = dt; 69 ((PC_Factor*)ilu)->info.dtcol = dtcol; 70 ((PC_Factor*)ilu)->info.dtcount = dtcount; 71 ((PC_Factor*)ilu)->info.fill = PETSC_DEFAULT; 72 PetscFunctionReturn(0); 73 } 74 EXTERN_C_END 75 76 EXTERN_C_BEGIN 77 #undef __FUNCT__ 78 #define __FUNCT__ "PCFactorSetReuseOrdering_ILU" 79 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering_ILU(PC pc,PetscTruth flag) 80 { 81 PC_ILU *ilu = (PC_ILU*)pc->data; 82 83 PetscFunctionBegin; 84 ilu->reuseordering = flag; 85 PetscFunctionReturn(0); 86 } 87 EXTERN_C_END 88 89 EXTERN_C_BEGIN 90 #undef __FUNCT__ 91 #define __FUNCT__ "PCFactorSetUseInPlace_ILU" 92 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseInPlace_ILU(PC pc) 93 { 94 PC_ILU *dir = (PC_ILU*)pc->data; 95 96 PetscFunctionBegin; 97 dir->inplace = PETSC_TRUE; 98 PetscFunctionReturn(0); 99 } 100 EXTERN_C_END 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "PCSetFromOptions_ILU" 104 static PetscErrorCode PCSetFromOptions_ILU(PC pc) 105 { 106 PetscErrorCode ierr; 107 PetscInt dtmax = 3,itmp; 108 PetscTruth flg,set; 109 PetscReal dt[3]; 110 char tname[256], solvertype[64]; 111 PC_ILU *ilu = (PC_ILU*)pc->data; 112 PetscFList ordlist; 113 PetscReal tol; 114 115 PetscFunctionBegin; 116 if (!MatOrderingRegisterAllCalled) {ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 117 ierr = PetscOptionsHead("ILU Options");CHKERRQ(ierr); 118 ierr = PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);CHKERRQ(ierr); 119 if (flg) ((PC_Factor*)ilu)->info.levels = itmp; 120 ierr = PetscOptionsTruth("-pc_factor_in_place","do factorization in place","PCFactorSetUseInPlace",ilu->inplace,&ilu->inplace,PETSC_NULL);CHKERRQ(ierr); 121 flg = PETSC_FALSE; 122 ierr = PetscOptionsTruth("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 123 ((PC_Factor*)ilu)->info.diagonal_fill = (double) flg; 124 ierr = PetscOptionsTruth("-pc_factor_reuse_fill","Reuse fill ratio from previous factorization","PCFactorSetReuseFill",ilu->reusefill,&ilu->reusefill,PETSC_NULL);CHKERRQ(ierr); 125 ierr = PetscOptionsTruth("-pc_factor_reuse_ordering","Reuse previous reordering","PCFactorSetReuseOrdering",ilu->reuseordering,&ilu->reuseordering,PETSC_NULL);CHKERRQ(ierr); 126 flg = PETSC_FALSE; 127 ierr = PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);CHKERRQ(ierr); 128 if (flg) { 129 ierr = PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);CHKERRQ(ierr); 130 } 131 ierr = PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",((PC_Factor*)ilu)->info.shiftnz,&((PC_Factor*)ilu)->info.shiftnz,0);CHKERRQ(ierr); 132 flg = (((PC_Factor*)ilu)->info.shiftpd > 0.0) ? PETSC_TRUE : PETSC_FALSE; 133 ierr = PetscOptionsTruth("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 134 ierr = PCFactorSetShiftPd(pc,flg);CHKERRQ(ierr); 135 ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)ilu)->info.zeropivot,&((PC_Factor*)ilu)->info.zeropivot,0);CHKERRQ(ierr); 136 137 dt[0] = ((PC_Factor*)ilu)->info.dt; 138 dt[1] = ((PC_Factor*)ilu)->info.dtcol; 139 dt[2] = ((PC_Factor*)ilu)->info.dtcount; 140 ierr = PetscOptionsRealArray("-pc_factor_use_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetUseDropTolerance",dt,&dtmax,&flg);CHKERRQ(ierr); 141 if (flg) { 142 ierr = PCFactorSetUseDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);CHKERRQ(ierr); 143 } 144 ierr = PetscOptionsReal("-pc_factor_fill","Expected fill in factorization","PCFactorSetFill",((PC_Factor*)ilu)->info.fill,&((PC_Factor*)ilu)->info.fill,&flg);CHKERRQ(ierr); 145 ierr = PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);CHKERRQ(ierr); 146 if (flg) { 147 tol = PETSC_DECIDE; 148 ierr = PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);CHKERRQ(ierr); 149 ierr = PCFactorReorderForNonzeroDiagonal(pc,tol);CHKERRQ(ierr); 150 } 151 152 ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr); 153 ierr = PetscOptionsList("-pc_factor_mat_ordering_type","Reorder to reduce nonzeros in ILU","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)ilu)->ordering,tname,256,&flg);CHKERRQ(ierr); 154 if (flg) { 155 ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr); 156 } 157 158 /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */ 159 ierr = PetscOptionsString("-pc_factor_mat_solver_package","Specific ILU solver to use","MatGetFactor",((PC_Factor*)ilu)->solvertype,solvertype,64,&flg);CHKERRQ(ierr); 160 if (flg) { 161 ierr = PCFactorSetMatSolverPackage(pc,solvertype);CHKERRQ(ierr); 162 } 163 164 flg = ((PC_Factor*)ilu)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE; 165 ierr = PetscOptionsTruth("-pc_factor_pivot_in_blocks","Pivot inside matrix blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",flg,&flg,&set);CHKERRQ(ierr); 166 if (set) { 167 ierr = PCFactorSetPivotInBlocks(pc,flg);CHKERRQ(ierr); 168 } 169 flg = PETSC_FALSE; 170 ierr = PetscOptionsTruth("-pc_factor_shift_in_blocks","Shift added to diagonal of block","PCFactorSetShiftInBlocks",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 171 if (flg) { 172 ierr = PCFactorSetShiftInBlocks(pc,(PetscReal)PETSC_DECIDE);CHKERRQ(ierr); 173 } 174 ierr = PetscOptionsReal("-pc_factor_shift_in_blocks","Shift added to diagonal of block","PCFactorSetShiftInBlocks",((PC_Factor*)ilu)->info.shiftinblocks,&((PC_Factor*)ilu)->info.shiftinblocks,0);CHKERRQ(ierr); 175 176 ierr = PetscOptionsTail();CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "PCView_ILU" 182 static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer) 183 { 184 PC_ILU *ilu = (PC_ILU*)pc->data; 185 PetscErrorCode ierr; 186 PetscTruth isstring,iascii; 187 188 PetscFunctionBegin; 189 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 190 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 191 if (iascii) { 192 if (ilu->usedt) { 193 ierr = PetscViewerASCIIPrintf(viewer," ILU: drop tolerance %G\n",((PC_Factor*)ilu)->info.dt);CHKERRQ(ierr); 194 ierr = PetscViewerASCIIPrintf(viewer," ILU: max nonzeros per row %D\n",(PetscInt)((PC_Factor*)ilu)->info.dtcount);CHKERRQ(ierr); 195 ierr = PetscViewerASCIIPrintf(viewer," ILU: column permutation tolerance %G\n",((PC_Factor*)ilu)->info.dtcol);CHKERRQ(ierr); 196 } else if (((PC_Factor*)ilu)->info.levels == 1) { 197 ierr = PetscViewerASCIIPrintf(viewer," ILU: %D level of fill\n",(PetscInt)((PC_Factor*)ilu)->info.levels);CHKERRQ(ierr); 198 } else { 199 ierr = PetscViewerASCIIPrintf(viewer," ILU: %D levels of fill\n",(PetscInt)((PC_Factor*)ilu)->info.levels);CHKERRQ(ierr); 200 } 201 ierr = PetscViewerASCIIPrintf(viewer," ILU: factor fill ratio allocated %G\n",((PC_Factor*)ilu)->info.fill);CHKERRQ(ierr); 202 ierr = PetscViewerASCIIPrintf(viewer," ILU: tolerance for zero pivot %G\n",((PC_Factor*)ilu)->info.zeropivot);CHKERRQ(ierr); 203 if (((PC_Factor*)ilu)->info.shiftpd) {ierr = PetscViewerASCIIPrintf(viewer," ILU: using Manteuffel shift\n");CHKERRQ(ierr);} 204 if (((PC_Factor*)ilu)->info.shiftnz) {ierr = PetscViewerASCIIPrintf(viewer," ILU: using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr);} 205 if (((PC_Factor*)ilu)->info.shiftinblocks) {ierr = PetscViewerASCIIPrintf(viewer," ILU: using diagonal shift on blocks to prevent zero pivot\n");CHKERRQ(ierr);} 206 if (ilu->inplace) {ierr = PetscViewerASCIIPrintf(viewer," in-place factorization\n");CHKERRQ(ierr);} 207 else {ierr = PetscViewerASCIIPrintf(viewer," out-of-place factorization\n");CHKERRQ(ierr);} 208 ierr = PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",((PC_Factor*)ilu)->ordering);CHKERRQ(ierr); 209 if (ilu->reusefill) {ierr = PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");CHKERRQ(ierr);} 210 if (ilu->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");CHKERRQ(ierr);} 211 if (((PC_Factor*)ilu)->fact) { 212 ierr = PetscViewerASCIIPrintf(viewer," ILU: factor fill ratio needed %G\n",ilu->actualfill);CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");CHKERRQ(ierr); 214 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 215 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 216 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 217 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 218 ierr = MatView(((PC_Factor*)ilu)->fact,viewer);CHKERRQ(ierr); 219 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 220 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 221 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 222 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 223 } 224 } else if (isstring) { 225 ierr = PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)((PC_Factor*)ilu)->info.levels,((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);CHKERRQ(ierr); 226 } else { 227 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCILU",((PetscObject)viewer)->type_name); 228 } 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "PCSetUp_ILU" 234 static PetscErrorCode PCSetUp_ILU(PC pc) 235 { 236 PetscErrorCode ierr; 237 PC_ILU *ilu = (PC_ILU*)pc->data; 238 MatInfo info; 239 240 PetscFunctionBegin; 241 if (ilu->inplace) { 242 CHKMEMQ; 243 if (!pc->setupcalled) { 244 245 /* In-place factorization only makes sense with the natural ordering, 246 so we only need to get the ordering once, even if nonzero structure changes */ 247 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr); 248 if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);} 249 if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);} 250 } 251 252 /* In place ILU only makes sense with fill factor of 1.0 because 253 cannot have levels of fill */ 254 ((PC_Factor*)ilu)->info.fill = 1.0; 255 ((PC_Factor*)ilu)->info.diagonal_fill = 0; 256 ierr = MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 257 CHKMEMQ; 258 ((PC_Factor*)ilu)->fact = pc->pmat; 259 } else if (ilu->usedt) { 260 if (!pc->setupcalled) { 261 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr); 262 CHKMEMQ; 263 ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ILUDT,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 264 if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);} 265 if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);} 266 ierr = MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 267 ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 268 } else if (pc->flag != SAME_NONZERO_PATTERN) { 269 CHKMEMQ; 270 ierr = MatDestroy(((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 271 CHKMEMQ; 272 if (!ilu->reuseordering) { 273 if (ilu->row) {ierr = ISDestroy(ilu->row);CHKERRQ(ierr);} 274 if (ilu->col) {ierr = ISDestroy(ilu->col);CHKERRQ(ierr);} 275 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr); 276 if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);} 277 if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);} 278 } 279 ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ILUDT,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 280 ierr = MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 281 ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 282 } else if (!ilu->reusefill) { 283 ierr = MatDestroy(((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 284 ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ILUDT,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 285 ierr = MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 286 ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 287 } else { 288 ierr = MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 289 } 290 } else { 291 if (!pc->setupcalled) { 292 /* first time in so compute reordering and symbolic factorization */ 293 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr); 294 if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);} 295 if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);} 296 /* Remove zeros along diagonal? */ 297 if (ilu->nonzerosalongdiagonal) { 298 ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr); 299 } 300 CHKMEMQ; 301 ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 302 CHKMEMQ; 303 ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 304 ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 305 ilu->actualfill = info.fill_ratio_needed; 306 ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 307 } else if (pc->flag != SAME_NONZERO_PATTERN) { 308 if (!ilu->reuseordering) { 309 /* compute a new ordering for the ILU */ 310 ierr = ISDestroy(ilu->row);CHKERRQ(ierr); 311 ierr = ISDestroy(ilu->col);CHKERRQ(ierr); 312 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr); 313 if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);} 314 if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);} 315 /* Remove zeros along diagonal? */ 316 if (ilu->nonzerosalongdiagonal) { 317 ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr); 318 } 319 } 320 ierr = MatDestroy(((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 321 ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 322 ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 323 ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 324 ilu->actualfill = info.fill_ratio_needed; 325 ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 326 } 327 CHKMEMQ; 328 ierr = MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 329 CHKMEMQ; 330 } 331 PetscFunctionReturn(0); 332 } 333 334 #undef __FUNCT__ 335 #define __FUNCT__ "PCDestroy_ILU" 336 static PetscErrorCode PCDestroy_ILU(PC pc) 337 { 338 PC_ILU *ilu = (PC_ILU*)pc->data; 339 PetscErrorCode ierr; 340 341 PetscFunctionBegin; 342 ierr = PCDestroy_ILU_Internal(pc);CHKERRQ(ierr); 343 ierr = PetscStrfree(((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr); 344 ierr = PetscStrfree(((PC_Factor*)ilu)->ordering);CHKERRQ(ierr); 345 ierr = PetscFree(ilu);CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348 349 #undef __FUNCT__ 350 #define __FUNCT__ "PCApply_ILU" 351 static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y) 352 { 353 PC_ILU *ilu = (PC_ILU*)pc->data; 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 ierr = MatSolve(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr); 358 PetscFunctionReturn(0); 359 } 360 361 #undef __FUNCT__ 362 #define __FUNCT__ "PCApplyTranspose_ILU" 363 static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y) 364 { 365 PC_ILU *ilu = (PC_ILU*)pc->data; 366 PetscErrorCode ierr; 367 368 PetscFunctionBegin; 369 ierr = MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr); 370 PetscFunctionReturn(0); 371 } 372 373 /*MC 374 PCILU - Incomplete factorization preconditioners. 375 376 Options Database Keys: 377 + -pc_factor_levels <k> - number of levels of fill for ILU(k) 378 . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for 379 its factorization (overwrites original matrix) 380 . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill 381 . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization 382 . -pc_factor_use_drop_tolerance <dt,dtcol,maxrowcount> - use Saad's drop tolerance ILUdt 383 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 384 . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal, 385 this decreases the chance of getting a zero pivot 386 . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 387 . -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger 388 than 1 the diagonal blocks are factored with partial pivoting (this increases the 389 stability of the ILU factorization 390 . -pc_factor_shift_in_blocks - adds a small diagonal to any block if it is singular during ILU factorization 391 . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default 392 - -pc_factor_shift_positive_definite true or false - Activate/Deactivate PCFactorSetShiftPd(); the value 393 is optional with true being the default 394 395 Level: beginner 396 397 Concepts: incomplete factorization 398 399 Notes: Only implemented for some matrix formats. (for parallel use you 400 must use MATMPIROWBS, see MatCreateMPIRowbs(), this supports only ILU(0) and this is not recommended 401 unless you really want a parallel ILU). 402 403 For BAIJ matrices this implements a point block ILU 404 405 References: 406 T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving 407 self-adjoint elliptic difference equations. SIAM J. Numer. Anal., 5:559--573, 1968. 408 409 T.A. Oliphant. An implicit numerical method for solving two-dimensional time-dependent dif- 410 fusion problems. Quart. Appl. Math., 19:221--229, 1961. 411 412 Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST 413 http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical 414 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 415 Science and Engineering, Kluwer, pp. 167--202. 416 417 418 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 419 PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), PCFactorSetUseDropTolerance(), 420 PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 421 PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(), 422 PCFactorSetShiftNonzero(),PCFactorSetShiftPd() 423 424 M*/ 425 426 EXTERN_C_BEGIN 427 #undef __FUNCT__ 428 #define __FUNCT__ "PCCreate_ILU" 429 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ILU(PC pc) 430 { 431 PetscErrorCode ierr; 432 PC_ILU *ilu; 433 434 PetscFunctionBegin; 435 ierr = PetscNewLog(pc,PC_ILU,&ilu);CHKERRQ(ierr); 436 437 ((PC_Factor*)ilu)->fact = 0; 438 ierr = MatFactorInfoInitialize(&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 439 ((PC_Factor*)ilu)->info.levels = 0; 440 ((PC_Factor*)ilu)->info.fill = 1.0; 441 ilu->col = 0; 442 ilu->row = 0; 443 ilu->inplace = PETSC_FALSE; 444 ierr = PetscStrallocpy(MAT_SOLVER_PETSC,&((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr); 445 ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)ilu)->ordering);CHKERRQ(ierr); 446 ilu->reuseordering = PETSC_FALSE; 447 ilu->usedt = PETSC_FALSE; 448 ((PC_Factor*)ilu)->info.dt = PETSC_DEFAULT; 449 ((PC_Factor*)ilu)->info.dtcount = PETSC_DEFAULT; 450 ((PC_Factor*)ilu)->info.dtcol = PETSC_DEFAULT; 451 ((PC_Factor*)ilu)->info.shiftnz = 1.e-12; 452 ((PC_Factor*)ilu)->info.shiftpd = 0.0; /* false */ 453 ((PC_Factor*)ilu)->info.zeropivot = 1.e-12; 454 ((PC_Factor*)ilu)->info.pivotinblocks = 1.0; 455 ((PC_Factor*)ilu)->info.shiftinblocks = 1.e-12; 456 ilu->reusefill = PETSC_FALSE; 457 ((PC_Factor*)ilu)->info.diagonal_fill = 0; 458 pc->data = (void*)ilu; 459 460 pc->ops->destroy = PCDestroy_ILU; 461 pc->ops->apply = PCApply_ILU; 462 pc->ops->applytranspose = PCApplyTranspose_ILU; 463 pc->ops->setup = PCSetUp_ILU; 464 pc->ops->setfromoptions = PCSetFromOptions_ILU; 465 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 466 pc->ops->view = PCView_ILU; 467 pc->ops->applyrichardson = 0; 468 469 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 470 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 471 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Factor", 472 PCFactorSetShiftNonzero_Factor);CHKERRQ(ierr); 473 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Factor", 474 PCFactorSetShiftPd_Factor);CHKERRQ(ierr); 475 476 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 477 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 478 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor", 479 PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 480 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseDropTolerance_C","PCFactorSetUseDropTolerance_ILU", 481 PCFactorSetUseDropTolerance_ILU);CHKERRQ(ierr); 482 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 483 PCFactorSetFill_Factor);CHKERRQ(ierr); 484 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 485 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 486 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_ILU", 487 PCFactorSetReuseOrdering_ILU);CHKERRQ(ierr); 488 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_ILU", 489 PCFactorSetReuseFill_ILU);CHKERRQ(ierr); 490 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor", 491 PCFactorSetLevels_Factor);CHKERRQ(ierr); 492 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_ILU", 493 PCFactorSetUseInPlace_ILU);CHKERRQ(ierr); 494 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C","PCFactorSetAllowDiagonalFill_Factor", 495 PCFactorSetAllowDiagonalFill_Factor);CHKERRQ(ierr); 496 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetPivotInBlocks_C","PCFactorSetPivotInBlocks_Factor", 497 PCFactorSetPivotInBlocks_Factor);CHKERRQ(ierr); 498 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftInBlocks_C","PCFactorSetShiftInBlocks_Factor", 499 PCFactorSetShiftInBlocks_Factor);CHKERRQ(ierr); 500 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C","PCFactorReorderForNonzeroDiagonal_ILU", 501 PCFactorReorderForNonzeroDiagonal_ILU);CHKERRQ(ierr); 502 PetscFunctionReturn(0); 503 } 504 EXTERN_C_END 505