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 } 102 103 /* In place ILU only makes sense with fill factor of 1.0 because 104 cannot have levels of fill */ 105 ((PC_Factor *)ilu)->info.fill = 1.0; 106 ((PC_Factor *)ilu)->info.diagonal_fill = 0.0; 107 108 PetscCall(MatILUFactor(pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info)); 109 PetscCall(MatFactorGetError(pc->pmat, &err)); 110 if (err) { /* Factor() fails */ 111 pc->failedreason = (PCFailedReason)err; 112 PetscFunctionReturn(0); 113 } 114 115 ((PC_Factor *)ilu)->fact = pc->pmat; 116 /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */ 117 PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &pc->matstate)); 118 } else { 119 if (!pc->setupcalled) { 120 /* first time in so compute reordering and symbolic factorization */ 121 PetscBool canuseordering; 122 if (!((PC_Factor *)ilu)->fact) { PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact)); } 123 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering)); 124 if (canuseordering) { 125 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 126 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col)); 127 /* Remove zeros along diagonal? */ 128 if (ilu->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col)); 129 } 130 PetscCall(MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info)); 131 PetscCall(MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info)); 132 ilu->hdr.actualfill = info.fill_ratio_needed; 133 } else if (pc->flag != SAME_NONZERO_PATTERN) { 134 if (!ilu->hdr.reuseordering) { 135 PetscBool canuseordering; 136 PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact)); 137 PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact)); 138 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering)); 139 if (canuseordering) { 140 /* compute a new ordering for the ILU */ 141 PetscCall(ISDestroy(&ilu->row)); 142 PetscCall(ISDestroy(&ilu->col)); 143 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 144 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col)); 145 /* Remove zeros along diagonal? */ 146 if (ilu->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col)); 147 } 148 } 149 PetscCall(MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info)); 150 PetscCall(MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info)); 151 ilu->hdr.actualfill = info.fill_ratio_needed; 152 } 153 PetscCall(MatFactorGetError(((PC_Factor *)ilu)->fact, &err)); 154 if (err) { /* FactorSymbolic() fails */ 155 pc->failedreason = (PCFailedReason)err; 156 PetscFunctionReturn(0); 157 } 158 159 PetscCall(MatLUFactorNumeric(((PC_Factor *)ilu)->fact, pc->pmat, &((PC_Factor *)ilu)->info)); 160 PetscCall(MatFactorGetError(((PC_Factor *)ilu)->fact, &err)); 161 if (err) { /* FactorNumeric() fails */ 162 pc->failedreason = (PCFailedReason)err; 163 } 164 } 165 166 PetscCall(PCFactorGetMatSolverType(pc, &stype)); 167 if (!stype) { 168 MatSolverType solverpackage; 169 PetscCall(MatFactorGetSolverType(((PC_Factor *)ilu)->fact, &solverpackage)); 170 PetscCall(PCFactorSetMatSolverType(pc, solverpackage)); 171 } 172 PetscFunctionReturn(0); 173 } 174 175 static PetscErrorCode PCDestroy_ILU(PC pc) { 176 PC_ILU *ilu = (PC_ILU *)pc->data; 177 178 PetscFunctionBegin; 179 PetscCall(PCReset_ILU(pc)); 180 PetscCall(PetscFree(((PC_Factor *)ilu)->solvertype)); 181 PetscCall(PetscFree(((PC_Factor *)ilu)->ordering)); 182 PetscCall(PetscFree(pc->data)); 183 PetscCall(PCFactorClearComposedFunctions(pc)); 184 PetscFunctionReturn(0); 185 } 186 187 static PetscErrorCode PCApply_ILU(PC pc, Vec x, Vec y) { 188 PC_ILU *ilu = (PC_ILU *)pc->data; 189 190 PetscFunctionBegin; 191 PetscCall(MatSolve(((PC_Factor *)ilu)->fact, x, y)); 192 PetscFunctionReturn(0); 193 } 194 195 static PetscErrorCode PCMatApply_ILU(PC pc, Mat X, Mat Y) { 196 PC_ILU *ilu = (PC_ILU *)pc->data; 197 198 PetscFunctionBegin; 199 PetscCall(MatMatSolve(((PC_Factor *)ilu)->fact, X, Y)); 200 PetscFunctionReturn(0); 201 } 202 203 static PetscErrorCode PCApplyTranspose_ILU(PC pc, Vec x, Vec y) { 204 PC_ILU *ilu = (PC_ILU *)pc->data; 205 206 PetscFunctionBegin; 207 PetscCall(MatSolveTranspose(((PC_Factor *)ilu)->fact, x, y)); 208 PetscFunctionReturn(0); 209 } 210 211 static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc, Vec x, Vec y) { 212 PC_ILU *icc = (PC_ILU *)pc->data; 213 214 PetscFunctionBegin; 215 PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y)); 216 PetscFunctionReturn(0); 217 } 218 219 static PetscErrorCode PCApplySymmetricRight_ILU(PC pc, Vec x, Vec y) { 220 PC_ILU *icc = (PC_ILU *)pc->data; 221 222 PetscFunctionBegin; 223 PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y)); 224 PetscFunctionReturn(0); 225 } 226 227 /*MC 228 PCILU - Incomplete factorization preconditioners. 229 230 Options Database Keys: 231 + -pc_factor_levels <k> - number of levels of fill for ILU(k) 232 . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for 233 its factorization (overwrites original matrix) 234 . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill 235 . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization 236 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 237 . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal, 238 this decreases the chance of getting a zero pivot 239 . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 240 - -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger 241 than 1 the diagonal blocks are factored with partial pivoting (this increases the 242 stability of the ILU factorization 243 244 Level: beginner 245 246 Notes: 247 Only implemented for some matrix format and sequential. For parallel see `PCHYPRE` for hypre's ILU 248 249 For `MATSEQBAIJ` matrices this implements a point block ILU 250 251 The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U) 252 even when the matrix is not symmetric since the U stores the diagonals of the factorization. 253 254 If you are using `MATSEQAIJCUSPARSE` matrices (or `MATMPIAIJCUSPARSE` matrices with block Jacobi), factorization 255 is never done on the GPU). 256 257 References: 258 + * - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving 259 self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968. 260 . * - T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961. 261 - * - TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS, 262 Chapter in Parallel Numerical 263 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 264 Science and Engineering, Kluwer. 265 266 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCLU`, `PCICC`, `PCCHOLESKY`, 267 `PCFactorSetZeroPivot()`, `PCFactorSetShiftSetType()`, `PCFactorSetAmount()`, 268 `PCFactorSetDropTolerance()`, `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`, 269 `PCFactorSetLevels()`, `PCFactorSetUseInPlace()`, `PCFactorSetAllowDiagonalFill()`, `PCFactorSetPivotInBlocks()`, 270 `PCFactorGetAllowDiagonalFill()`, `PCFactorGetUseInPlace()` 271 M*/ 272 273 PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc) { 274 PC_ILU *ilu; 275 276 PetscFunctionBegin; 277 PetscCall(PetscNew(&ilu)); 278 pc->data = (void *)ilu; 279 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ILU)); 280 281 ((PC_Factor *)ilu)->info.levels = 0.; 282 ((PC_Factor *)ilu)->info.fill = 1.0; 283 ilu->col = NULL; 284 ilu->row = NULL; 285 ((PC_Factor *)ilu)->info.dt = PETSC_DEFAULT; 286 ((PC_Factor *)ilu)->info.dtcount = PETSC_DEFAULT; 287 ((PC_Factor *)ilu)->info.dtcol = PETSC_DEFAULT; 288 289 pc->ops->reset = PCReset_ILU; 290 pc->ops->destroy = PCDestroy_ILU; 291 pc->ops->apply = PCApply_ILU; 292 pc->ops->matapply = PCMatApply_ILU; 293 pc->ops->applytranspose = PCApplyTranspose_ILU; 294 pc->ops->setup = PCSetUp_ILU; 295 pc->ops->setfromoptions = PCSetFromOptions_ILU; 296 pc->ops->view = PCView_Factor; 297 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU; 298 pc->ops->applysymmetricright = PCApplySymmetricRight_ILU; 299 pc->ops->applyrichardson = NULL; 300 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetDropTolerance_C", PCFactorSetDropTolerance_ILU)); 301 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_ILU)); 302 PetscFunctionReturn(0); 303 } 304