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 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 PCApplyTranspose_ILU(PC pc,Vec x,Vec y) 207 { 208 PC_ILU *ilu = (PC_ILU*)pc->data; 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 ierr = MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y) 217 { 218 PetscErrorCode ierr; 219 PC_ILU *icc = (PC_ILU*)pc->data; 220 221 PetscFunctionBegin; 222 ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 223 PetscFunctionReturn(0); 224 } 225 226 static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y) 227 { 228 PetscErrorCode ierr; 229 PC_ILU *icc = (PC_ILU*)pc->data; 230 231 PetscFunctionBegin; 232 ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 233 PetscFunctionReturn(0); 234 } 235 236 /*MC 237 PCILU - Incomplete factorization preconditioners. 238 239 Options Database Keys: 240 + -pc_factor_levels <k> - number of levels of fill for ILU(k) 241 . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for 242 its factorization (overwrites original matrix) 243 . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill 244 . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization 245 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 246 . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal, 247 this decreases the chance of getting a zero pivot 248 . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 249 - -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger 250 than 1 the diagonal blocks are factored with partial pivoting (this increases the 251 stability of the ILU factorization 252 253 Level: beginner 254 255 Notes: 256 Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU) 257 258 For BAIJ matrices this implements a point block ILU 259 260 The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U) 261 even when the matrix is not symmetric since the U stores the diagonals of the factorization. 262 263 If you are using MATSEQAIJCUSPARSE matrices (or MATMPIAIJCUSPARESE matrices with block Jacobi), factorization 264 is never done on the GPU). 265 266 References: 267 + 1. - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving 268 self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968. 269 . 2. - T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961. 270 - 3. - TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS, 271 Chapter in Parallel Numerical 272 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 273 Science and Engineering, Kluwer. 274 275 276 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 277 PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(), 278 PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 279 PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(), 280 PCFactorGetAllowDiagonalFill(), PCFactorGetUseInPlace() 281 282 M*/ 283 284 PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc) 285 { 286 PetscErrorCode ierr; 287 PC_ILU *ilu; 288 289 PetscFunctionBegin; 290 ierr = PetscNewLog(pc,&ilu);CHKERRQ(ierr); 291 pc->data = (void*)ilu; 292 ierr = PCFactorInitialize(pc);CHKERRQ(ierr); 293 294 ((PC_Factor*)ilu)->factortype = MAT_FACTOR_ILU; 295 ((PC_Factor*)ilu)->info.levels = 0.; 296 ((PC_Factor*)ilu)->info.fill = 1.0; 297 ilu->col = 0; 298 ilu->row = 0; 299 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)ilu)->ordering);CHKERRQ(ierr); 300 ((PC_Factor*)ilu)->info.dt = PETSC_DEFAULT; 301 ((PC_Factor*)ilu)->info.dtcount = PETSC_DEFAULT; 302 ((PC_Factor*)ilu)->info.dtcol = PETSC_DEFAULT; 303 304 pc->ops->reset = PCReset_ILU; 305 pc->ops->destroy = PCDestroy_ILU; 306 pc->ops->apply = PCApply_ILU; 307 pc->ops->applytranspose = PCApplyTranspose_ILU; 308 pc->ops->setup = PCSetUp_ILU; 309 pc->ops->setfromoptions = PCSetFromOptions_ILU; 310 pc->ops->view = PCView_Factor; 311 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU; 312 pc->ops->applysymmetricright = PCApplySymmetricRight_ILU; 313 pc->ops->applyrichardson = 0; 314 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);CHKERRQ(ierr); 315 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_ILU);CHKERRQ(ierr); 316 PetscFunctionReturn(0); 317 } 318