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. 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 PCFactorSetShiftPd(pc,PETSC_FALSE); 166 to turn off the shift. 167 168 References: 169 Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST 170 http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical 171 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 172 Science and Engineering, Kluwer, pp. 167--202. 173 174 175 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 176 PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(),PCFactorSetShiftInBlocks(), 177 PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 178 PCFactorSetLevels() 179 180 M*/ 181 182 EXTERN_C_BEGIN 183 #undef __FUNCT__ 184 #define __FUNCT__ "PCCreate_ICC" 185 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ICC(PC pc) 186 { 187 PetscErrorCode ierr; 188 PC_ICC *icc; 189 190 PetscFunctionBegin; 191 ierr = PetscNewLog(pc,PC_ICC,&icc);CHKERRQ(ierr); 192 193 ((PC_Factor*)icc)->fact = 0; 194 ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 195 ierr = MatFactorInfoInitialize(&((PC_Factor*)icc)->info);CHKERRQ(ierr); 196 ((PC_Factor*)icc)->info.levels = 0; 197 ((PC_Factor*)icc)->info.fill = 1.0; 198 icc->implctx = 0; 199 200 ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT; 201 ((PC_Factor*)icc)->info.shiftnz = 0.0; 202 ((PC_Factor*)icc)->info.shiftpd = 1.0; /* true */ 203 ((PC_Factor*)icc)->info.zeropivot = 1.e-12; 204 pc->data = (void*)icc; 205 206 pc->ops->apply = PCApply_ICC; 207 pc->ops->setup = PCSetup_ICC; 208 pc->ops->destroy = PCDestroy_ICC; 209 pc->ops->setfromoptions = PCSetFromOptions_ICC; 210 pc->ops->view = PCView_ICC; 211 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 212 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 213 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 214 215 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 216 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 217 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 218 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 219 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Factor", 220 PCFactorSetShiftNonzero_Factor);CHKERRQ(ierr); 221 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Factor", 222 PCFactorSetShiftPd_Factor);CHKERRQ(ierr); 223 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftInBlocks_C","PCFactorSetShiftInBlocks_Factor", 224 PCFactorSetShiftInBlocks_Factor);CHKERRQ(ierr); 225 226 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor", 227 PCFactorSetLevels_Factor);CHKERRQ(ierr); 228 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 229 PCFactorSetFill_Factor);CHKERRQ(ierr); 230 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 231 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 EXTERN_C_END 235 236 237