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 const MatSolverPackage stype; 13 MatFactorError err; 14 15 PetscFunctionBegin; 16 ierr = MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);CHKERRQ(ierr); 17 18 ierr = MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);CHKERRQ(ierr); 19 if (!pc->setupcalled) { 20 if (!((PC_Factor*)icc)->fact) { 21 ierr = MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);CHKERRQ(ierr); 22 } 23 ierr = MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 24 } else if (pc->flag != SAME_NONZERO_PATTERN) { 25 ierr = MatDestroy(&((PC_Factor*)icc)->fact);CHKERRQ(ierr); 26 ierr = MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);CHKERRQ(ierr); 27 ierr = MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 28 } 29 ierr = MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 30 icc->hdr.actualfill = info.fill_ratio_needed; 31 32 ierr = ISDestroy(&cperm);CHKERRQ(ierr); 33 ierr = ISDestroy(&perm);CHKERRQ(ierr); 34 35 ierr = MatFactorGetError(((PC_Factor*)icc)->fact,&err);CHKERRQ(ierr); 36 if (err) { /* FactorSymbolic() fails */ 37 pc->failedreason = (PCFailedReason)err; 38 PetscFunctionReturn(0); 39 } 40 41 ierr = MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 42 ierr = MatFactorGetError(((PC_Factor*)icc)->fact,&err);CHKERRQ(ierr); 43 if (err) { /* FactorNumeric() fails */ 44 pc->failedreason = (PCFailedReason)err; 45 } 46 47 ierr = PCFactorGetMatSolverPackage(pc,&stype);CHKERRQ(ierr); 48 if (!stype) { 49 const MatSolverPackage solverpackage; 50 ierr = MatFactorGetSolverPackage(((PC_Factor*)icc)->fact,&solverpackage);CHKERRQ(ierr); 51 ierr = PCFactorSetMatSolverPackage(pc,solverpackage);CHKERRQ(ierr); 52 } 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "PCReset_ICC" 58 static PetscErrorCode PCReset_ICC(PC pc) 59 { 60 PC_ICC *icc = (PC_ICC*)pc->data; 61 PetscErrorCode ierr; 62 63 PetscFunctionBegin; 64 ierr = MatDestroy(&((PC_Factor*)icc)->fact);CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "PCDestroy_ICC" 70 static PetscErrorCode PCDestroy_ICC(PC pc) 71 { 72 PC_ICC *icc = (PC_ICC*)pc->data; 73 PetscErrorCode ierr; 74 75 PetscFunctionBegin; 76 ierr = PCReset_ICC(pc);CHKERRQ(ierr); 77 ierr = PetscFree(((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 78 ierr = PetscFree(((PC_Factor*)icc)->solvertype);CHKERRQ(ierr); 79 ierr = PetscFree(pc->data);CHKERRQ(ierr); 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "PCApply_ICC" 85 static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y) 86 { 87 PC_ICC *icc = (PC_ICC*)pc->data; 88 PetscErrorCode ierr; 89 90 PetscFunctionBegin; 91 ierr = MatSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "PCApplySymmetricLeft_ICC" 97 static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y) 98 { 99 PetscErrorCode ierr; 100 PC_ICC *icc = (PC_ICC*)pc->data; 101 102 PetscFunctionBegin; 103 ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 104 PetscFunctionReturn(0); 105 } 106 107 #undef __FUNCT__ 108 #define __FUNCT__ "PCApplySymmetricRight_ICC" 109 static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y) 110 { 111 PetscErrorCode ierr; 112 PC_ICC *icc = (PC_ICC*)pc->data; 113 114 PetscFunctionBegin; 115 ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "PCSetFromOptions_ICC" 121 static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc) 122 { 123 PC_ICC *icc = (PC_ICC*)pc->data; 124 PetscBool flg; 125 PetscErrorCode ierr; 126 /* PetscReal dt[3];*/ 127 128 PetscFunctionBegin; 129 ierr = PetscOptionsHead(PetscOptionsObject,"ICC Options");CHKERRQ(ierr); 130 ierr = PCSetFromOptions_Factor(PetscOptionsObject,pc);CHKERRQ(ierr); 131 132 ierr = PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);CHKERRQ(ierr); 133 /*dt[0] = ((PC_Factor*)icc)->info.dt; 134 dt[1] = ((PC_Factor*)icc)->info.dtcol; 135 dt[2] = ((PC_Factor*)icc)->info.dtcount; 136 PetscInt dtmax = 3; 137 ierr = PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);CHKERRQ(ierr); 138 if (flg) { 139 ierr = PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);CHKERRQ(ierr); 140 } 141 */ 142 ierr = PetscOptionsTail();CHKERRQ(ierr); 143 PetscFunctionReturn(0); 144 } 145 146 #undef __FUNCT__ 147 #define __FUNCT__ "PCView_ICC" 148 static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer) 149 { 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 ierr = PCView_Factor(pc,viewer);CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt); 158 159 /*MC 160 PCICC - Incomplete Cholesky factorization preconditioners. 161 162 Options Database Keys: 163 + -pc_factor_levels <k> - number of levels of fill for ICC(k) 164 . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for 165 its factorization (overwrites original matrix) 166 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 167 - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 168 169 Level: beginner 170 171 Concepts: incomplete Cholesky factorization 172 173 Notes: Only implemented for some matrix formats. Not implemented in parallel. 174 175 For BAIJ matrices this implements a point block ICC. 176 177 The Manteuffel shift is only implemented for matrices with block size 1 178 179 By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE); 180 to turn off the shift. 181 182 References: 183 . 1. - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, 184 Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 185 Science and Engineering, Kluwer. 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 pc->data = (void*)icc; 205 ierr = PCFactorInitialize(pc);CHKERRQ(ierr); 206 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 207 208 ((PC_Factor*)icc)->factortype = MAT_FACTOR_ICC; 209 ((PC_Factor*)icc)->info.fill = 1.0; 210 ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT; 211 ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE; 212 213 pc->ops->apply = PCApply_ICC; 214 pc->ops->applytranspose = PCApply_ICC; 215 pc->ops->setup = PCSetUp_ICC; 216 pc->ops->reset = PCReset_ICC; 217 pc->ops->destroy = PCDestroy_ICC; 218 pc->ops->setfromoptions = PCSetFromOptions_ICC; 219 pc->ops->view = PCView_ICC; 220 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 221 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 222 PetscFunctionReturn(0); 223 } 224 225 226