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