1 2 /* 3 Defines a ILU factorization preconditioner for any Mat implementation 4 */ 5 #include <../src/ksp/pc/impls/factor/ilu/ilu.h> /*I "petscpc.h" I*/ 6 7 PetscErrorCode PCFactorReorderForNonzeroDiagonal_ILU(PC pc, PetscReal z) { 8 PC_ILU *ilu = (PC_ILU *)pc->data; 9 10 PetscFunctionBegin; 11 ilu->nonzerosalongdiagonal = PETSC_TRUE; 12 if (z == PETSC_DECIDE) ilu->nonzerosalongdiagonaltol = 1.e-10; 13 else ilu->nonzerosalongdiagonaltol = z; 14 PetscFunctionReturn(0); 15 } 16 17 PetscErrorCode PCReset_ILU(PC pc) { 18 PC_ILU *ilu = (PC_ILU *)pc->data; 19 20 PetscFunctionBegin; 21 if (!ilu->hdr.inplace) PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact)); 22 if (ilu->row && ilu->col && ilu->row != ilu->col) PetscCall(ISDestroy(&ilu->row)); 23 PetscCall(ISDestroy(&ilu->col)); 24 PetscFunctionReturn(0); 25 } 26 27 PetscErrorCode PCFactorSetDropTolerance_ILU(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount) { 28 PC_ILU *ilu = (PC_ILU *)pc->data; 29 30 PetscFunctionBegin; 31 if (pc->setupcalled && (((PC_Factor *)ilu)->info.dt != dt || ((PC_Factor *)ilu)->info.dtcol != dtcol || ((PC_Factor *)ilu)->info.dtcount != dtcount)) { 32 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot change drop tolerance after using PC"); 33 } 34 ((PC_Factor *)ilu)->info.dt = dt; 35 ((PC_Factor *)ilu)->info.dtcol = dtcol; 36 ((PC_Factor *)ilu)->info.dtcount = dtcount; 37 ((PC_Factor *)ilu)->info.usedt = 1.0; 38 PetscFunctionReturn(0); 39 } 40 41 static PetscErrorCode PCSetFromOptions_ILU(PC pc, PetscOptionItems *PetscOptionsObject) { 42 PetscInt itmp; 43 PetscBool flg, set; 44 PC_ILU *ilu = (PC_ILU *)pc->data; 45 PetscReal tol; 46 47 PetscFunctionBegin; 48 PetscOptionsHeadBegin(PetscOptionsObject, "ILU Options"); 49 PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject)); 50 51 PetscCall(PetscOptionsInt("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", (PetscInt)((PC_Factor *)ilu)->info.levels, &itmp, &flg)); 52 if (flg) ((PC_Factor *)ilu)->info.levels = itmp; 53 54 PetscCall(PetscOptionsBool("-pc_factor_diagonal_fill", "Allow fill into empty diagonal entry", "PCFactorSetAllowDiagonalFill", ((PC_Factor *)ilu)->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE, &flg, &set)); 55 if (set) ((PC_Factor *)ilu)->info.diagonal_fill = (PetscReal)flg; 56 PetscCall(PetscOptionsName("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", &flg)); 57 if (flg) { 58 tol = PETSC_DECIDE; 59 PetscCall(PetscOptionsReal("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", ilu->nonzerosalongdiagonaltol, &tol, NULL)); 60 PetscCall(PCFactorReorderForNonzeroDiagonal(pc, tol)); 61 } 62 63 PetscOptionsHeadEnd(); 64 PetscFunctionReturn(0); 65 } 66 67 static PetscErrorCode PCSetUp_ILU(PC pc) { 68 PC_ILU *ilu = (PC_ILU *)pc->data; 69 MatInfo info; 70 PetscBool flg; 71 MatSolverType stype; 72 MatFactorError err; 73 const char *prefix; 74 75 PetscFunctionBegin; 76 pc->failedreason = PC_NOERROR; 77 /* ugly hack to change default, since it is not support by some matrix types */ 78 if (((PC_Factor *)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) { 79 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg)); 80 if (!flg) { 81 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPIAIJ, &flg)); 82 if (!flg) { 83 ((PC_Factor *)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS; 84 PetscCall(PetscInfo(pc, "Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO\n")); 85 } 86 } 87 } 88 89 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 90 PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix)); 91 92 PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure)); 93 if (ilu->hdr.inplace) { 94 if (!pc->setupcalled) { 95 /* In-place factorization only makes sense with the natural ordering, 96 so we only need to get the ordering once, even if nonzero structure changes */ 97 /* Should not get the ordering if the factorization routine does not use it, but do not yet have access to the factor matrix */ 98 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 99 PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact)); 100 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col)); 101 if (ilu->row) PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)ilu->row)); 102 if (ilu->col) PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)ilu->col)); 103 } 104 105 /* In place ILU only makes sense with fill factor of 1.0 because 106 cannot have levels of fill */ 107 ((PC_Factor *)ilu)->info.fill = 1.0; 108 ((PC_Factor *)ilu)->info.diagonal_fill = 0.0; 109 110 PetscCall(MatILUFactor(pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info)); 111 PetscCall(MatFactorGetError(pc->pmat, &err)); 112 if (err) { /* Factor() fails */ 113 pc->failedreason = (PCFailedReason)err; 114 PetscFunctionReturn(0); 115 } 116 117 ((PC_Factor *)ilu)->fact = pc->pmat; 118 /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */ 119 PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &pc->matstate)); 120 } else { 121 if (!pc->setupcalled) { 122 /* first time in so compute reordering and symbolic factorization */ 123 PetscBool canuseordering; 124 if (!((PC_Factor *)ilu)->fact) { 125 PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact)); 126 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)((PC_Factor *)ilu)->fact)); 127 } 128 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering)); 129 if (canuseordering) { 130 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 131 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col)); 132 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)ilu->row)); 133 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)ilu->col)); 134 /* Remove zeros along diagonal? */ 135 if (ilu->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col)); 136 } 137 PetscCall(MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info)); 138 PetscCall(MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info)); 139 ilu->hdr.actualfill = info.fill_ratio_needed; 140 } else if (pc->flag != SAME_NONZERO_PATTERN) { 141 if (!ilu->hdr.reuseordering) { 142 PetscBool canuseordering; 143 PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact)); 144 PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact)); 145 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)((PC_Factor *)ilu)->fact)); 146 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering)); 147 if (canuseordering) { 148 /* compute a new ordering for the ILU */ 149 PetscCall(ISDestroy(&ilu->row)); 150 PetscCall(ISDestroy(&ilu->col)); 151 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 152 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col)); 153 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)ilu->row)); 154 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)ilu->col)); 155 /* Remove zeros along diagonal? */ 156 if (ilu->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col)); 157 } 158 } 159 PetscCall(MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info)); 160 PetscCall(MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info)); 161 ilu->hdr.actualfill = info.fill_ratio_needed; 162 } 163 PetscCall(MatFactorGetError(((PC_Factor *)ilu)->fact, &err)); 164 if (err) { /* FactorSymbolic() fails */ 165 pc->failedreason = (PCFailedReason)err; 166 PetscFunctionReturn(0); 167 } 168 169 PetscCall(MatLUFactorNumeric(((PC_Factor *)ilu)->fact, pc->pmat, &((PC_Factor *)ilu)->info)); 170 PetscCall(MatFactorGetError(((PC_Factor *)ilu)->fact, &err)); 171 if (err) { /* FactorNumeric() fails */ 172 pc->failedreason = (PCFailedReason)err; 173 } 174 } 175 176 PetscCall(PCFactorGetMatSolverType(pc, &stype)); 177 if (!stype) { 178 MatSolverType solverpackage; 179 PetscCall(MatFactorGetSolverType(((PC_Factor *)ilu)->fact, &solverpackage)); 180 PetscCall(PCFactorSetMatSolverType(pc, solverpackage)); 181 } 182 PetscFunctionReturn(0); 183 } 184 185 static PetscErrorCode PCDestroy_ILU(PC pc) { 186 PC_ILU *ilu = (PC_ILU *)pc->data; 187 188 PetscFunctionBegin; 189 PetscCall(PCReset_ILU(pc)); 190 PetscCall(PetscFree(((PC_Factor *)ilu)->solvertype)); 191 PetscCall(PetscFree(((PC_Factor *)ilu)->ordering)); 192 PetscCall(PetscFree(pc->data)); 193 PetscCall(PCFactorClearComposedFunctions(pc)); 194 PetscFunctionReturn(0); 195 } 196 197 static PetscErrorCode PCApply_ILU(PC pc, Vec x, Vec y) { 198 PC_ILU *ilu = (PC_ILU *)pc->data; 199 200 PetscFunctionBegin; 201 PetscCall(MatSolve(((PC_Factor *)ilu)->fact, x, y)); 202 PetscFunctionReturn(0); 203 } 204 205 static PetscErrorCode PCMatApply_ILU(PC pc, Mat X, Mat Y) { 206 PC_ILU *ilu = (PC_ILU *)pc->data; 207 208 PetscFunctionBegin; 209 PetscCall(MatMatSolve(((PC_Factor *)ilu)->fact, X, Y)); 210 PetscFunctionReturn(0); 211 } 212 213 static PetscErrorCode PCApplyTranspose_ILU(PC pc, Vec x, Vec y) { 214 PC_ILU *ilu = (PC_ILU *)pc->data; 215 216 PetscFunctionBegin; 217 PetscCall(MatSolveTranspose(((PC_Factor *)ilu)->fact, x, y)); 218 PetscFunctionReturn(0); 219 } 220 221 static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc, Vec x, Vec y) { 222 PC_ILU *icc = (PC_ILU *)pc->data; 223 224 PetscFunctionBegin; 225 PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y)); 226 PetscFunctionReturn(0); 227 } 228 229 static PetscErrorCode PCApplySymmetricRight_ILU(PC pc, Vec x, Vec y) { 230 PC_ILU *icc = (PC_ILU *)pc->data; 231 232 PetscFunctionBegin; 233 PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y)); 234 PetscFunctionReturn(0); 235 } 236 237 /*MC 238 PCILU - Incomplete factorization preconditioners. 239 240 Options Database Keys: 241 + -pc_factor_levels <k> - number of levels of fill for ILU(k) 242 . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for 243 its factorization (overwrites original matrix) 244 . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill 245 . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization 246 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 247 . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal, 248 this decreases the chance of getting a zero pivot 249 . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 250 - -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger 251 than 1 the diagonal blocks are factored with partial pivoting (this increases the 252 stability of the ILU factorization 253 254 Level: beginner 255 256 Notes: 257 Only implemented for some matrix format and sequential. For parallel see `PCHYPRE` for hypre's ILU 258 259 For `MATSEQBAIJ` matrices this implements a point block ILU 260 261 The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U) 262 even when the matrix is not symmetric since the U stores the diagonals of the factorization. 263 264 If you are using `MATSEQAIJCUSPARSE` matrices (or `MATMPIAIJCUSPARSE` matrices with block Jacobi), factorization 265 is never done on the GPU). 266 267 References: 268 + * - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving 269 self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968. 270 . * - T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961. 271 - * - TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS, 272 Chapter in Parallel Numerical 273 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 274 Science and Engineering, Kluwer. 275 276 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCLU`, `PCICC`, `PCCHOLESKY`, 277 `PCFactorSetZeroPivot()`, `PCFactorSetShiftSetType()`, `PCFactorSetAmount()`, 278 `PCFactorSetDropTolerance()`, `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`, 279 `PCFactorSetLevels()`, `PCFactorSetUseInPlace()`, `PCFactorSetAllowDiagonalFill()`, `PCFactorSetPivotInBlocks()`, 280 `PCFactorGetAllowDiagonalFill()`, `PCFactorGetUseInPlace()` 281 M*/ 282 283 PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc) { 284 PC_ILU *ilu; 285 286 PetscFunctionBegin; 287 PetscCall(PetscNewLog(pc, &ilu)); 288 pc->data = (void *)ilu; 289 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ILU)); 290 291 ((PC_Factor *)ilu)->info.levels = 0.; 292 ((PC_Factor *)ilu)->info.fill = 1.0; 293 ilu->col = NULL; 294 ilu->row = NULL; 295 ((PC_Factor *)ilu)->info.dt = PETSC_DEFAULT; 296 ((PC_Factor *)ilu)->info.dtcount = PETSC_DEFAULT; 297 ((PC_Factor *)ilu)->info.dtcol = PETSC_DEFAULT; 298 299 pc->ops->reset = PCReset_ILU; 300 pc->ops->destroy = PCDestroy_ILU; 301 pc->ops->apply = PCApply_ILU; 302 pc->ops->matapply = PCMatApply_ILU; 303 pc->ops->applytranspose = PCApplyTranspose_ILU; 304 pc->ops->setup = PCSetUp_ILU; 305 pc->ops->setfromoptions = PCSetFromOptions_ILU; 306 pc->ops->view = PCView_Factor; 307 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU; 308 pc->ops->applysymmetricright = PCApplySymmetricRight_ILU; 309 pc->ops->applyrichardson = NULL; 310 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetDropTolerance_C", PCFactorSetDropTolerance_ILU)); 311 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_ILU)); 312 PetscFunctionReturn(0); 313 } 314