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