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 ierr = MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,& ((PC_Factor*)icc)->fact);CHKERRQ(ierr); 19 ierr = MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 20 } else if (pc->flag != SAME_NONZERO_PATTERN) { 21 ierr = MatDestroy(((PC_Factor*)icc)->fact);CHKERRQ(ierr); 22 ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);CHKERRQ(ierr); 23 ierr = MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 24 } 25 ierr = MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 26 icc->actualfill = info.fill_ratio_needed; 27 28 ierr = ISDestroy(cperm);CHKERRQ(ierr); 29 ierr = ISDestroy(perm);CHKERRQ(ierr); 30 ierr = MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 #undef __FUNCT__ 35 #define __FUNCT__ "PCDestroy_ICC" 36 static PetscErrorCode PCDestroy_ICC(PC pc) 37 { 38 PC_ICC *icc = (PC_ICC*)pc->data; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 if (((PC_Factor*)icc)->fact) {ierr = MatDestroy(((PC_Factor*)icc)->fact);CHKERRQ(ierr);} 43 ierr = PetscStrfree(((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 44 ierr = PetscStrfree(((PC_Factor*)icc)->solvertype);CHKERRQ(ierr); 45 ierr = PetscFree(icc);CHKERRQ(ierr); 46 PetscFunctionReturn(0); 47 } 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "PCApply_ICC" 51 static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y) 52 { 53 PC_ICC *icc = (PC_ICC*)pc->data; 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 ierr = MatSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "PCApplySymmetricLeft_ICC" 63 static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y) 64 { 65 PetscErrorCode ierr; 66 PC_ICC *icc = (PC_ICC*)pc->data; 67 68 PetscFunctionBegin; 69 ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "PCApplySymmetricRight_ICC" 75 static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y) 76 { 77 PetscErrorCode ierr; 78 PC_ICC *icc = (PC_ICC*)pc->data; 79 80 PetscFunctionBegin; 81 ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 #undef __FUNCT__ 86 #define __FUNCT__ "PCSetFromOptions_ICC" 87 static PetscErrorCode PCSetFromOptions_ICC(PC pc) 88 { 89 PC_ICC *icc = (PC_ICC*)pc->data; 90 PetscTruth flg; 91 PetscInt dtmax = 3; 92 PetscReal dt[3]; 93 PetscErrorCode ierr; 94 95 PetscFunctionBegin; 96 ierr = PetscOptionsHead("ICC Options");CHKERRQ(ierr); 97 ierr = PCSetFromOptions_Factor(pc);CHKERRQ(ierr); 98 99 ierr = PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);CHKERRQ(ierr); 100 dt[0] = ((PC_Factor*)icc)->info.dt; 101 dt[1] = ((PC_Factor*)icc)->info.dtcol; 102 dt[2] = ((PC_Factor*)icc)->info.dtcount; 103 ierr = PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);CHKERRQ(ierr); 104 if (flg) { 105 ierr = PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);CHKERRQ(ierr); 106 } 107 108 ierr = PetscOptionsTail();CHKERRQ(ierr); 109 PetscFunctionReturn(0); 110 } 111 112 #undef __FUNCT__ 113 #define __FUNCT__ "PCView_ICC" 114 static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer) 115 { 116 PC_ICC *icc = (PC_ICC*)pc->data; 117 PetscErrorCode ierr; 118 PetscTruth isstring,iascii; 119 120 PetscFunctionBegin; 121 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 122 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 123 if (iascii) { 124 if (((PC_Factor*)icc)->info.levels == 1) { 125 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D level of fill\n",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr); 126 } else { 127 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D levels of fill\n",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr); 128 } 129 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio allocated %G, ordering used %s\n",((PC_Factor*)icc)->info.fill,((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 130 131 if (((PC_Factor*)icc)->info.shifttype==MAT_SHIFT_POSITIVE_DEFINITE) { 132 ierr = PetscViewerASCIIPrintf(viewer," ICC: using Manteuffel shift\n");CHKERRQ(ierr); 133 } 134 if (((PC_Factor*)icc)->info.shifttype==MAT_SHIFT_NONZERO) { 135 ierr = PetscViewerASCIIPrintf(viewer," ICC: using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr); 136 } 137 if (((PC_Factor*)icc)->info.shifttype==MAT_SHIFT_INBLOCKS) { 138 ierr = PetscViewerASCIIPrintf(viewer," ICC: using diagonal shift on blocks to prevent zero pivot\n");CHKERRQ(ierr); 139 } 140 141 if (((PC_Factor*)icc)->fact) { 142 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio needed %G\n",icc->actualfill);CHKERRQ(ierr); 143 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");CHKERRQ(ierr); 144 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 145 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 146 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 147 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 148 ierr = MatView(((PC_Factor*)icc)->fact,viewer);CHKERRQ(ierr); 149 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 150 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 151 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 152 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 153 } 154 } else if (isstring) { 155 ierr = PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr);CHKERRQ(ierr); 156 } else { 157 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name); 158 } 159 PetscFunctionReturn(0); 160 } 161 162 EXTERN_C_BEGIN 163 extern PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt); 164 EXTERN_C_END 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 Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST 191 http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical 192 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 193 Science and Engineering, Kluwer, pp. 167--202. 194 195 196 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 197 PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(), 198 PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 199 PCFactorSetLevels() 200 201 M*/ 202 203 EXTERN_C_BEGIN 204 #undef __FUNCT__ 205 #define __FUNCT__ "PCCreate_ICC" 206 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ICC(PC pc) 207 { 208 PetscErrorCode ierr; 209 PC_ICC *icc; 210 211 PetscFunctionBegin; 212 ierr = PetscNewLog(pc,PC_ICC,&icc);CHKERRQ(ierr); 213 214 ((PC_Factor*)icc)->fact = 0; 215 ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 216 ierr = PetscStrallocpy(MAT_SOLVER_PETSC,&((PC_Factor*)icc)->solvertype);CHKERRQ(ierr); 217 ierr = MatFactorInfoInitialize(&((PC_Factor*)icc)->info);CHKERRQ(ierr); 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 = MAT_SHIFT_POSITIVE_DEFINITE; 224 ((PC_Factor*)icc)->info.shiftamount = 0.0; 225 ((PC_Factor*)icc)->info.zeropivot = 1.e-12; 226 227 pc->data = (void*)icc; 228 pc->ops->apply = PCApply_ICC; 229 pc->ops->setup = PCSetup_ICC; 230 pc->ops->destroy = PCDestroy_ICC; 231 pc->ops->setfromoptions = PCSetFromOptions_ICC; 232 pc->ops->view = PCView_ICC; 233 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 234 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 235 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 236 237 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 238 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 239 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 240 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 241 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor", 242 PCFactorSetShiftType_Factor);CHKERRQ(ierr); 243 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor", 244 PCFactorSetShiftAmount_Factor);CHKERRQ(ierr); 245 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor", 246 PCFactorSetLevels_Factor);CHKERRQ(ierr); 247 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 248 PCFactorSetFill_Factor);CHKERRQ(ierr); 249 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 250 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 251 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor", 252 PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 253 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetDropTolerance_C","PCFactorSetDropTolerance_ILU", 254 PCFactorSetDropTolerance_ILU);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 EXTERN_C_END 258 259 260