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 formats. (for parallel see PCHYPRE for hypre's ILU) 258 259 For BAIJ 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`, 277 `PCFactorSetZeroPivot()`, `PCFactorSetShiftSetType()`, `PCFactorSetAmount()`, 278 `PCFactorSetDropTolerance()`, `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`, 279 `PCFactorSetLevels()`, `PCFactorSetUseInPlace()`, `PCFactorSetAllowDiagonalFill()`, `PCFactorSetPivotInBlocks()`, 280 `PCFactorGetAllowDiagonalFill()`, `PCFactorGetUseInPlace()` 281 282 M*/ 283 284 PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc) { 285 PC_ILU *ilu; 286 287 PetscFunctionBegin; 288 PetscCall(PetscNewLog(pc, &ilu)); 289 pc->data = (void *)ilu; 290 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ILU)); 291 292 ((PC_Factor *)ilu)->info.levels = 0.; 293 ((PC_Factor *)ilu)->info.fill = 1.0; 294 ilu->col = NULL; 295 ilu->row = NULL; 296 ((PC_Factor *)ilu)->info.dt = PETSC_DEFAULT; 297 ((PC_Factor *)ilu)->info.dtcount = PETSC_DEFAULT; 298 ((PC_Factor *)ilu)->info.dtcol = PETSC_DEFAULT; 299 300 pc->ops->reset = PCReset_ILU; 301 pc->ops->destroy = PCDestroy_ILU; 302 pc->ops->apply = PCApply_ILU; 303 pc->ops->matapply = PCMatApply_ILU; 304 pc->ops->applytranspose = PCApplyTranspose_ILU; 305 pc->ops->setup = PCSetUp_ILU; 306 pc->ops->setfromoptions = PCSetFromOptions_ILU; 307 pc->ops->view = PCView_Factor; 308 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU; 309 pc->ops->applysymmetricright = PCApplySymmetricRight_ILU; 310 pc->ops->applyrichardson = NULL; 311 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetDropTolerance_C", PCFactorSetDropTolerance_ILU)); 312 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_ILU)); 313 PetscFunctionReturn(0); 314 } 315