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