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 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr); 103 if (ilu->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);CHKERRQ(ierr);} 104 if (ilu->col) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);CHKERRQ(ierr);} 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 ierr = MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 113 ierr = MatFactorGetError(pc->pmat,&err);CHKERRQ(ierr); 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 ierr = PetscObjectStateGet((PetscObject)pc->pmat,&pc->matstate);CHKERRQ(ierr); 122 } else { 123 if (!pc->setupcalled) { 124 /* first time in so compute reordering and symbolic factorization */ 125 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr); 126 if (ilu->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);CHKERRQ(ierr);} 127 if (ilu->col) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);CHKERRQ(ierr);} 128 /* Remove zeros along diagonal? */ 129 if (ilu->nonzerosalongdiagonal) { 130 ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr); 131 } 132 if (!((PC_Factor*)ilu)->fact) { 133 ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 134 } 135 ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 136 ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 137 ilu->hdr.actualfill = info.fill_ratio_needed; 138 139 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 140 } else if (pc->flag != SAME_NONZERO_PATTERN) { 141 if (!ilu->hdr.reuseordering) { 142 /* compute a new ordering for the ILU */ 143 ierr = ISDestroy(&ilu->row);CHKERRQ(ierr); 144 ierr = ISDestroy(&ilu->col);CHKERRQ(ierr); 145 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr); 146 if (ilu->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);CHKERRQ(ierr);} 147 if (ilu->col) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);CHKERRQ(ierr);} 148 /* Remove zeros along diagonal? */ 149 if (ilu->nonzerosalongdiagonal) { 150 ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr); 151 } 152 } 153 ierr = MatDestroy(&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 154 ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 155 ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 156 ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 157 ilu->hdr.actualfill = info.fill_ratio_needed; 158 159 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);CHKERRQ(ierr); 160 } 161 ierr = MatFactorGetError(((PC_Factor*)ilu)->fact,&err);CHKERRQ(ierr); 162 if (err) { /* FactorSymbolic() fails */ 163 pc->failedreason = (PCFailedReason)err; 164 PetscFunctionReturn(0); 165 } 166 167 ierr = MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);CHKERRQ(ierr); 168 ierr = MatFactorGetError(((PC_Factor*)ilu)->fact,&err);CHKERRQ(ierr); 169 if (err) { /* FactorNumeric() fails */ 170 pc->failedreason = (PCFailedReason)err; 171 } 172 } 173 174 ierr = PCFactorGetMatSolverType(pc,&stype);CHKERRQ(ierr); 175 if (!stype) { 176 MatSolverType solverpackage; 177 ierr = MatFactorGetSolverType(((PC_Factor*)ilu)->fact,&solverpackage);CHKERRQ(ierr); 178 ierr = PCFactorSetMatSolverType(pc,solverpackage);CHKERRQ(ierr); 179 } 180 PetscFunctionReturn(0); 181 } 182 183 static PetscErrorCode PCDestroy_ILU(PC pc) 184 { 185 PC_ILU *ilu = (PC_ILU*)pc->data; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = PCReset_ILU(pc);CHKERRQ(ierr); 190 ierr = PetscFree(((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr); 191 ierr = PetscFree(((PC_Factor*)ilu)->ordering);CHKERRQ(ierr); 192 ierr = PetscFree(pc->data);CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y) 197 { 198 PC_ILU *ilu = (PC_ILU*)pc->data; 199 PetscErrorCode ierr; 200 201 PetscFunctionBegin; 202 ierr = MatSolve(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 static PetscErrorCode PCMatApply_ILU(PC pc,Mat X,Mat Y) 207 { 208 PC_ILU *ilu = (PC_ILU*)pc->data; 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 ierr = MatMatSolve(((PC_Factor*)ilu)->fact,X,Y);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y) 217 { 218 PC_ILU *ilu = (PC_ILU*)pc->data; 219 PetscErrorCode ierr; 220 221 PetscFunctionBegin; 222 ierr = MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr); 223 PetscFunctionReturn(0); 224 } 225 226 static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y) 227 { 228 PetscErrorCode ierr; 229 PC_ILU *icc = (PC_ILU*)pc->data; 230 231 PetscFunctionBegin; 232 ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 233 PetscFunctionReturn(0); 234 } 235 236 static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y) 237 { 238 PetscErrorCode ierr; 239 PC_ILU *icc = (PC_ILU*)pc->data; 240 241 PetscFunctionBegin; 242 ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 /*MC 247 PCILU - Incomplete factorization preconditioners. 248 249 Options Database Keys: 250 + -pc_factor_levels <k> - number of levels of fill for ILU(k) 251 . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for 252 its factorization (overwrites original matrix) 253 . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill 254 . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization 255 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 256 . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal, 257 this decreases the chance of getting a zero pivot 258 . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 259 - -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger 260 than 1 the diagonal blocks are factored with partial pivoting (this increases the 261 stability of the ILU factorization 262 263 Level: beginner 264 265 Notes: 266 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 MATMPIAIJCUSPARSE 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 = NULL; 308 ilu->row = NULL; 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->matapply = PCMatApply_ILU; 318 pc->ops->applytranspose = PCApplyTranspose_ILU; 319 pc->ops->setup = PCSetUp_ILU; 320 pc->ops->setfromoptions = PCSetFromOptions_ILU; 321 pc->ops->view = PCView_Factor; 322 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU; 323 pc->ops->applysymmetricright = PCApplySymmetricRight_ILU; 324 pc->ops->applyrichardson = NULL; 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