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