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