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