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