1 2 #include <../src/ksp/pc/impls/factor/icc/icc.h> /*I "petscpc.h" I*/ 3 4 static PetscErrorCode PCSetUp_ICC(PC pc) { 5 PC_ICC *icc = (PC_ICC *)pc->data; 6 IS perm = NULL, cperm = NULL; 7 MatInfo info; 8 MatSolverType stype; 9 MatFactorError err; 10 PetscBool canuseordering; 11 const char *prefix; 12 13 PetscFunctionBegin; 14 pc->failedreason = PC_NOERROR; 15 16 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 17 PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix)); 18 19 PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure)); 20 if (!pc->setupcalled) { 21 if (!((PC_Factor *)icc)->fact) { PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, MAT_FACTOR_ICC, &((PC_Factor *)icc)->fact)); } 22 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)icc)->fact, &canuseordering)); 23 if (canuseordering) { 24 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 25 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)icc)->ordering, &perm, &cperm)); 26 } 27 PetscCall(MatICCFactorSymbolic(((PC_Factor *)icc)->fact, pc->pmat, perm, &((PC_Factor *)icc)->info)); 28 } else if (pc->flag != SAME_NONZERO_PATTERN) { 29 PetscBool canuseordering; 30 PetscCall(MatDestroy(&((PC_Factor *)icc)->fact)); 31 PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, MAT_FACTOR_ICC, &((PC_Factor *)icc)->fact)); 32 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)icc)->fact, &canuseordering)); 33 if (canuseordering) { 34 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 35 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)icc)->ordering, &perm, &cperm)); 36 } 37 PetscCall(MatICCFactorSymbolic(((PC_Factor *)icc)->fact, pc->pmat, perm, &((PC_Factor *)icc)->info)); 38 } 39 PetscCall(MatGetInfo(((PC_Factor *)icc)->fact, MAT_LOCAL, &info)); 40 icc->hdr.actualfill = info.fill_ratio_needed; 41 42 PetscCall(ISDestroy(&cperm)); 43 PetscCall(ISDestroy(&perm)); 44 45 PetscCall(MatFactorGetError(((PC_Factor *)icc)->fact, &err)); 46 if (err) { /* FactorSymbolic() fails */ 47 pc->failedreason = (PCFailedReason)err; 48 PetscFunctionReturn(0); 49 } 50 51 PetscCall(MatCholeskyFactorNumeric(((PC_Factor *)icc)->fact, pc->pmat, &((PC_Factor *)icc)->info)); 52 PetscCall(MatFactorGetError(((PC_Factor *)icc)->fact, &err)); 53 if (err) { /* FactorNumeric() fails */ 54 pc->failedreason = (PCFailedReason)err; 55 } 56 57 PetscCall(PCFactorGetMatSolverType(pc, &stype)); 58 if (!stype) { 59 MatSolverType solverpackage; 60 PetscCall(MatFactorGetSolverType(((PC_Factor *)icc)->fact, &solverpackage)); 61 PetscCall(PCFactorSetMatSolverType(pc, solverpackage)); 62 } 63 PetscFunctionReturn(0); 64 } 65 66 static PetscErrorCode PCReset_ICC(PC pc) { 67 PC_ICC *icc = (PC_ICC *)pc->data; 68 69 PetscFunctionBegin; 70 PetscCall(MatDestroy(&((PC_Factor *)icc)->fact)); 71 PetscFunctionReturn(0); 72 } 73 74 static PetscErrorCode PCDestroy_ICC(PC pc) { 75 PC_ICC *icc = (PC_ICC *)pc->data; 76 77 PetscFunctionBegin; 78 PetscCall(PCReset_ICC(pc)); 79 PetscCall(PetscFree(((PC_Factor *)icc)->ordering)); 80 PetscCall(PetscFree(((PC_Factor *)icc)->solvertype)); 81 PetscCall(PCFactorClearComposedFunctions(pc)); 82 PetscCall(PetscFree(pc->data)); 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode PCApply_ICC(PC pc, Vec x, Vec y) { 87 PC_ICC *icc = (PC_ICC *)pc->data; 88 89 PetscFunctionBegin; 90 PetscCall(MatSolve(((PC_Factor *)icc)->fact, x, y)); 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode PCMatApply_ICC(PC pc, Mat X, Mat Y) { 95 PC_ICC *icc = (PC_ICC *)pc->data; 96 97 PetscFunctionBegin; 98 PetscCall(MatMatSolve(((PC_Factor *)icc)->fact, X, Y)); 99 PetscFunctionReturn(0); 100 } 101 102 static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc, Vec x, Vec y) { 103 PC_ICC *icc = (PC_ICC *)pc->data; 104 105 PetscFunctionBegin; 106 PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y)); 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode PCApplySymmetricRight_ICC(PC pc, Vec x, Vec y) { 111 PC_ICC *icc = (PC_ICC *)pc->data; 112 113 PetscFunctionBegin; 114 PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y)); 115 PetscFunctionReturn(0); 116 } 117 118 static PetscErrorCode PCSetFromOptions_ICC(PC pc, PetscOptionItems *PetscOptionsObject) { 119 PC_ICC *icc = (PC_ICC *)pc->data; 120 PetscBool flg; 121 /* PetscReal dt[3];*/ 122 123 PetscFunctionBegin; 124 PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options"); 125 PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject)); 126 127 PetscCall(PetscOptionsReal("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", ((PC_Factor *)icc)->info.levels, &((PC_Factor *)icc)->info.levels, &flg)); 128 /*dt[0] = ((PC_Factor*)icc)->info.dt; 129 dt[1] = ((PC_Factor*)icc)->info.dtcol; 130 dt[2] = ((PC_Factor*)icc)->info.dtcount; 131 PetscInt dtmax = 3; 132 PetscCall(PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg)); 133 if (flg) { 134 PetscCall(PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2])); 135 } 136 */ 137 PetscOptionsHeadEnd(); 138 PetscFunctionReturn(0); 139 } 140 141 extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC, PetscReal, PetscReal, PetscInt); 142 143 /*MC 144 PCICC - Incomplete Cholesky factorization preconditioners. 145 146 Options Database Keys: 147 + -pc_factor_levels <k> - number of levels of fill for ICC(k) 148 . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for 149 its factorization (overwrites original matrix) 150 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 151 - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 152 153 Level: beginner 154 155 Notes: 156 Only implemented for some matrix formats. Not implemented in parallel. 157 158 For BAIJ matrices this implements a point block ICC. 159 160 The Manteuffel shift is only implemented for matrices with block size 1 161 162 By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE); 163 to turn off the shift. 164 165 References: 166 . * - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, 167 Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 168 Science and Engineering, Kluwer. 169 170 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, 171 `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`, 172 `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`, 173 `PCFactorSetLevels()` 174 175 M*/ 176 177 PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc) { 178 PC_ICC *icc; 179 180 PetscFunctionBegin; 181 PetscCall(PetscNewLog(pc, &icc)); 182 pc->data = (void *)icc; 183 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ICC)); 184 185 ((PC_Factor *)icc)->info.fill = 1.0; 186 ((PC_Factor *)icc)->info.dtcol = PETSC_DEFAULT; 187 ((PC_Factor *)icc)->info.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE; 188 189 pc->ops->apply = PCApply_ICC; 190 pc->ops->matapply = PCMatApply_ICC; 191 pc->ops->applytranspose = PCApply_ICC; 192 pc->ops->setup = PCSetUp_ICC; 193 pc->ops->reset = PCReset_ICC; 194 pc->ops->destroy = PCDestroy_ICC; 195 pc->ops->setfromoptions = PCSetFromOptions_ICC; 196 pc->ops->view = PCView_Factor; 197 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 198 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 199 PetscFunctionReturn(0); 200 } 201