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,MAT_SOLVER_PETSC,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 = PetscFree(icc);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "PCApply_ICC" 50 static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y) 51 { 52 PC_ICC *icc = (PC_ICC*)pc->data; 53 PetscErrorCode ierr; 54 55 PetscFunctionBegin; 56 ierr = MatSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 57 PetscFunctionReturn(0); 58 } 59 60 #undef __FUNCT__ 61 #define __FUNCT__ "PCApplySymmetricLeft_ICC" 62 static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y) 63 { 64 PetscErrorCode ierr; 65 PC_ICC *icc = (PC_ICC*)pc->data; 66 67 PetscFunctionBegin; 68 ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "PCApplySymmetricRight_ICC" 74 static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y) 75 { 76 PetscErrorCode ierr; 77 PC_ICC *icc = (PC_ICC*)pc->data; 78 79 PetscFunctionBegin; 80 ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "PCSetFromOptions_ICC" 86 static PetscErrorCode PCSetFromOptions_ICC(PC pc) 87 { 88 PC_ICC *icc = (PC_ICC*)pc->data; 89 PetscTruth flg; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 ierr = PetscOptionsHead("ICC Options");CHKERRQ(ierr); 94 ierr = PCSetFromOptions_Factor(pc);CHKERRQ(ierr); 95 96 ierr = PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);CHKERRQ(ierr); 97 ierr = PetscOptionsTail();CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "PCView_ICC" 103 static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer) 104 { 105 PC_ICC *icc = (PC_ICC*)pc->data; 106 PetscErrorCode ierr; 107 PetscTruth isstring,iascii; 108 109 PetscFunctionBegin; 110 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 111 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 112 if (iascii) { 113 if (((PC_Factor*)icc)->info.levels == 1) { 114 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D level of fill\n",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr); 115 } else { 116 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D levels of fill\n",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr); 117 } 118 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio allocated %G, ordering used %s\n",((PC_Factor*)icc)->info.fill,((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 119 if (((PC_Factor*)icc)->info.shiftpd) {ierr = PetscViewerASCIIPrintf(viewer," ICC: using Manteuffel shift\n");CHKERRQ(ierr);} 120 if (((PC_Factor*)icc)->info.shiftnz) {ierr = PetscViewerASCIIPrintf(viewer," ICC: using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr);} 121 if (((PC_Factor*)icc)->fact) { 122 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio needed %G\n",icc->actualfill);CHKERRQ(ierr); 123 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");CHKERRQ(ierr); 124 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 125 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 126 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 127 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 128 ierr = MatView(((PC_Factor*)icc)->fact,viewer);CHKERRQ(ierr); 129 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 130 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 131 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 132 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 133 } 134 } else if (isstring) { 135 ierr = PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr);CHKERRQ(ierr); 136 } else { 137 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name); 138 } 139 PetscFunctionReturn(0); 140 } 141 142 /*MC 143 PCICC - Incomplete Cholesky factorization preconditioners. 144 145 Options Database Keys: 146 + -pc_factor_levels <k> - number of levels of fill for ICC(k) 147 . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for 148 its factorization (overwrites original matrix) 149 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 150 . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 151 . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default 152 - -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value 153 is optional with PETSC_TRUE being the default 154 155 Level: beginner 156 157 Concepts: incomplete Cholesky factorization 158 159 Notes: Only implemented for some matrix formats. Not implemented in parallel (for parallel use you 160 must use MATMPIROWBS, see MatCreateMPIRowbs(), this supports only ICC(0) and this is not recommended 161 unless you really want a parallel ICC). 162 163 For BAIJ matrices this implements a point block ICC. 164 165 The Manteuffel shift is only implemented for matrices with block size 1 166 167 By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftPd(pc,PETSC_FALSE); 168 to turn off the shift. 169 170 References: 171 Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST 172 http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical 173 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 174 Science and Engineering, Kluwer, pp. 167--202. 175 176 177 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 178 PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(),PCFactorSetShiftInBlocks(), 179 PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 180 PCFactorSetLevels() 181 182 M*/ 183 184 EXTERN_C_BEGIN 185 #undef __FUNCT__ 186 #define __FUNCT__ "PCCreate_ICC" 187 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ICC(PC pc) 188 { 189 PetscErrorCode ierr; 190 PC_ICC *icc; 191 192 PetscFunctionBegin; 193 ierr = PetscNewLog(pc,PC_ICC,&icc);CHKERRQ(ierr); 194 195 ((PC_Factor*)icc)->fact = 0; 196 ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 197 ierr = MatFactorInfoInitialize(&((PC_Factor*)icc)->info);CHKERRQ(ierr); 198 ((PC_Factor*)icc)->info.levels = 0; 199 ((PC_Factor*)icc)->info.fill = 1.0; 200 icc->implctx = 0; 201 202 ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT; 203 ((PC_Factor*)icc)->info.shiftnz = 0.0; 204 ((PC_Factor*)icc)->info.shiftpd = 1.0; /* true */ 205 ((PC_Factor*)icc)->info.zeropivot = 1.e-12; 206 pc->data = (void*)icc; 207 208 pc->ops->apply = PCApply_ICC; 209 pc->ops->setup = PCSetup_ICC; 210 pc->ops->destroy = PCDestroy_ICC; 211 pc->ops->setfromoptions = PCSetFromOptions_ICC; 212 pc->ops->view = PCView_ICC; 213 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 214 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 215 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 216 217 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 218 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 219 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 220 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 221 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Factor", 222 PCFactorSetShiftNonzero_Factor);CHKERRQ(ierr); 223 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Factor", 224 PCFactorSetShiftPd_Factor);CHKERRQ(ierr); 225 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftInBlocks_C","PCFactorSetShiftInBlocks_Factor", 226 PCFactorSetShiftInBlocks_Factor);CHKERRQ(ierr); 227 228 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor", 229 PCFactorSetLevels_Factor);CHKERRQ(ierr); 230 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 231 PCFactorSetFill_Factor);CHKERRQ(ierr); 232 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 233 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 234 PetscFunctionReturn(0); 235 } 236 EXTERN_C_END 237 238 239