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 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: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU) 267 268 For BAIJ matrices this implements a point block ILU 269 270 The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U) 271 even when the matrix is not symmetric since the U stores the diagonals of the factorization. 272 273 If you are using MATSEQAIJCUSPARSE matrices (or MATMPIAIJCUSPARESE matrices with block Jacobi), factorization 274 is never done on the GPU). 275 276 References: 277 + 1. - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving 278 self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968. 279 . 2. - T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961. 280 - 3. - TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS, 281 Chapter in Parallel Numerical 282 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 283 Science and Engineering, Kluwer. 284 285 286 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 287 PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(), 288 PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 289 PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(), 290 PCFactorGetAllowDiagonalFill(), PCFactorGetUseInPlace() 291 292 M*/ 293 294 PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc) 295 { 296 PetscErrorCode ierr; 297 PC_ILU *ilu; 298 299 PetscFunctionBegin; 300 ierr = PetscNewLog(pc,&ilu);CHKERRQ(ierr); 301 pc->data = (void*)ilu; 302 ierr = PCFactorInitialize(pc);CHKERRQ(ierr); 303 304 ((PC_Factor*)ilu)->factortype = MAT_FACTOR_ILU; 305 ((PC_Factor*)ilu)->info.levels = 0.; 306 ((PC_Factor*)ilu)->info.fill = 1.0; 307 ilu->col = 0; 308 ilu->row = 0; 309 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)ilu)->ordering);CHKERRQ(ierr); 310 ((PC_Factor*)ilu)->info.dt = PETSC_DEFAULT; 311 ((PC_Factor*)ilu)->info.dtcount = PETSC_DEFAULT; 312 ((PC_Factor*)ilu)->info.dtcol = PETSC_DEFAULT; 313 314 pc->ops->reset = PCReset_ILU; 315 pc->ops->destroy = PCDestroy_ILU; 316 pc->ops->apply = PCApply_ILU; 317 pc->ops->applytranspose = PCApplyTranspose_ILU; 318 pc->ops->setup = PCSetUp_ILU; 319 pc->ops->setfromoptions = PCSetFromOptions_ILU; 320 pc->ops->view = PCView_ILU; 321 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU; 322 pc->ops->applysymmetricright = PCApplySymmetricRight_ILU; 323 pc->ops->applyrichardson = 0; 324 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);CHKERRQ(ierr); 325 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_ILU);CHKERRQ(ierr); 326 PetscFunctionReturn(0); 327 } 328