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