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