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