1 2 #include <../src/ksp/pc/impls/factor/icc/icc.h> /*I "petscpc.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "PCSetUp_ICC" 6 static PetscErrorCode PCSetUp_ICC(PC pc) 7 { 8 PC_ICC *icc = (PC_ICC*)pc->data; 9 IS perm,cperm; 10 PetscErrorCode ierr; 11 MatInfo info; 12 Mat F; 13 14 PetscFunctionBegin; 15 ierr = MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);CHKERRQ(ierr); 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 = MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 23 } else if (pc->flag != SAME_NONZERO_PATTERN) { 24 ierr = MatDestroy(&((PC_Factor*)icc)->fact);CHKERRQ(ierr); 25 ierr = MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);CHKERRQ(ierr); 26 ierr = MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 27 } 28 ierr = MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 29 icc->actualfill = info.fill_ratio_needed; 30 31 ierr = ISDestroy(&cperm);CHKERRQ(ierr); 32 ierr = ISDestroy(&perm);CHKERRQ(ierr); 33 34 F = ((PC_Factor*)icc)->fact; 35 if (F->errortype) { /* FactorSymbolic() fails */ 36 pc->failedreason = (PCFailedReason)F->errortype; 37 PetscFunctionReturn(0); 38 } 39 40 ierr = MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 41 if (F->errortype) { /* FactorNumeric() fails */ 42 pc->failedreason = (PCFailedReason)F->errortype; 43 } 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "PCReset_ICC" 49 static PetscErrorCode PCReset_ICC(PC pc) 50 { 51 PC_ICC *icc = (PC_ICC*)pc->data; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 ierr = MatDestroy(&((PC_Factor*)icc)->fact);CHKERRQ(ierr); 56 PetscFunctionReturn(0); 57 } 58 59 #undef __FUNCT__ 60 #define __FUNCT__ "PCDestroy_ICC" 61 static PetscErrorCode PCDestroy_ICC(PC pc) 62 { 63 PC_ICC *icc = (PC_ICC*)pc->data; 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 ierr = PCReset_ICC(pc);CHKERRQ(ierr); 68 ierr = PetscFree(((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 69 ierr = PetscFree(((PC_Factor*)icc)->solvertype);CHKERRQ(ierr); 70 ierr = PetscFree(pc->data);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "PCApply_ICC" 76 static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y) 77 { 78 PC_ICC *icc = (PC_ICC*)pc->data; 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 ierr = MatSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 #undef __FUNCT__ 87 #define __FUNCT__ "PCApplySymmetricLeft_ICC" 88 static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y) 89 { 90 PetscErrorCode ierr; 91 PC_ICC *icc = (PC_ICC*)pc->data; 92 93 PetscFunctionBegin; 94 ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "PCApplySymmetricRight_ICC" 100 static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y) 101 { 102 PetscErrorCode ierr; 103 PC_ICC *icc = (PC_ICC*)pc->data; 104 105 PetscFunctionBegin; 106 ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "PCSetFromOptions_ICC" 112 static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc) 113 { 114 PC_ICC *icc = (PC_ICC*)pc->data; 115 PetscBool flg; 116 PetscErrorCode ierr; 117 /* PetscReal dt[3];*/ 118 119 PetscFunctionBegin; 120 ierr = PetscOptionsHead(PetscOptionsObject,"ICC Options");CHKERRQ(ierr); 121 ierr = PCSetFromOptions_Factor(PetscOptionsObject,pc);CHKERRQ(ierr); 122 123 ierr = PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);CHKERRQ(ierr); 124 /*dt[0] = ((PC_Factor*)icc)->info.dt; 125 dt[1] = ((PC_Factor*)icc)->info.dtcol; 126 dt[2] = ((PC_Factor*)icc)->info.dtcount; 127 PetscInt dtmax = 3; 128 ierr = PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);CHKERRQ(ierr); 129 if (flg) { 130 ierr = PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);CHKERRQ(ierr); 131 } 132 */ 133 ierr = PetscOptionsTail();CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 #undef __FUNCT__ 138 #define __FUNCT__ "PCView_ICC" 139 static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer) 140 { 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 ierr = PCView_Factor(pc,viewer);CHKERRQ(ierr); 145 PetscFunctionReturn(0); 146 } 147 148 extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt); 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "PCFactorGetUseInPlace_ICC" 152 PetscErrorCode PCFactorGetUseInPlace_ICC(PC pc,PetscBool *flg) 153 { 154 PetscFunctionBegin; 155 *flg = PETSC_FALSE; 156 PetscFunctionReturn(0); 157 } 158 159 /*MC 160 PCICC - Incomplete Cholesky factorization preconditioners. 161 162 Options Database Keys: 163 + -pc_factor_levels <k> - number of levels of fill for ICC(k) 164 . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for 165 its factorization (overwrites original matrix) 166 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 167 - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 168 169 Level: beginner 170 171 Concepts: incomplete Cholesky factorization 172 173 Notes: Only implemented for some matrix formats. Not implemented in parallel. 174 175 For BAIJ matrices this implements a point block ICC. 176 177 The Manteuffel shift is only implemented for matrices with block size 1 178 179 By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE); 180 to turn off the shift. 181 182 References: 183 TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, 184 http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical 185 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 186 Science and Engineering, Kluwer, pp. 167--202. 187 188 189 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 190 PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(), 191 PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 192 PCFactorSetLevels() 193 194 M*/ 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "PCCreate_ICC" 198 PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc) 199 { 200 PetscErrorCode ierr; 201 PC_ICC *icc; 202 203 PetscFunctionBegin; 204 ierr = PetscNewLog(pc,&icc);CHKERRQ(ierr); 205 206 ((PC_Factor*)icc)->fact = 0; 207 208 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 209 ierr = PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)icc)->solvertype);CHKERRQ(ierr); 210 ierr = MatFactorInfoInitialize(&((PC_Factor*)icc)->info);CHKERRQ(ierr); 211 212 ((PC_Factor*)icc)->factortype = MAT_FACTOR_ICC; 213 ((PC_Factor*)icc)->info.levels = 0.; 214 ((PC_Factor*)icc)->info.fill = 1.0; 215 icc->implctx = 0; 216 217 ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT; 218 ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE; 219 ((PC_Factor*)icc)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON; 220 ((PC_Factor*)icc)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON; 221 222 pc->data = (void*)icc; 223 pc->ops->apply = PCApply_ICC; 224 pc->ops->applytranspose = PCApply_ICC; 225 pc->ops->setup = PCSetUp_ICC; 226 pc->ops->reset = PCReset_ICC; 227 pc->ops->destroy = PCDestroy_ICC; 228 pc->ops->setfromoptions = PCSetFromOptions_ICC; 229 pc->ops->view = PCView_ICC; 230 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 231 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 232 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 233 234 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr); 235 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 236 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 237 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);CHKERRQ(ierr); 238 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);CHKERRQ(ierr); 239 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C",PCFactorSetLevels_Factor);CHKERRQ(ierr); 240 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetLevels_C",PCFactorGetLevels_Factor);CHKERRQ(ierr); 241 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);CHKERRQ(ierr); 242 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 243 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 244 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);CHKERRQ(ierr); 245 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_ICC);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 249 250