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