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