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