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