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