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