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