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 char tname[256]; 90 PetscTruth flg; 91 PetscErrorCode ierr; 92 PetscFList ordlist; 93 94 PetscFunctionBegin; 95 if (!MatOrderingRegisterAllCalled) {ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 96 ierr = PetscOptionsHead("ICC Options");CHKERRQ(ierr); 97 ierr = PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);CHKERRQ(ierr); 98 ierr = PetscOptionsReal("-pc_factor_fill","Expected fill in factorization","PCFactorSetFill",((PC_Factor*)icc)->info.fill,&((PC_Factor*)icc)->info.fill,&flg);CHKERRQ(ierr); 99 ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr); 100 ierr = PetscOptionsList("-pc_factor_mat_ordering_type","Reorder to reduce nonzeros in ICC","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)icc)->ordering,tname,256,&flg);CHKERRQ(ierr); 101 if (flg) { 102 ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr); 103 } 104 flg = PETSC_FALSE; 105 ierr = PetscOptionsTruth("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 106 if (flg) { 107 ierr = PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);CHKERRQ(ierr); 108 } 109 ierr = PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",((PC_Factor*)icc)->info.shiftnz,&((PC_Factor*)icc)->info.shiftnz,0);CHKERRQ(ierr); 110 flg = (((PC_Factor*)icc)->info.shiftpd > 0.0) ? PETSC_TRUE : PETSC_FALSE; 111 ierr = PetscOptionsTruth("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 112 ierr = PCFactorSetShiftPd(pc,flg);CHKERRQ(ierr); 113 ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)icc)->info.zeropivot,&((PC_Factor*)icc)->info.zeropivot,0);CHKERRQ(ierr); 114 115 ierr = PetscOptionsTail();CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "PCView_ICC" 121 static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer) 122 { 123 PC_ICC *icc = (PC_ICC*)pc->data; 124 PetscErrorCode ierr; 125 PetscTruth isstring,iascii; 126 127 PetscFunctionBegin; 128 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 129 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 130 if (iascii) { 131 if (((PC_Factor*)icc)->info.levels == 1) { 132 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D level of fill\n",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr); 133 } else { 134 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D levels of fill\n",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr); 135 } 136 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio allocated %G, ordering used %s\n",((PC_Factor*)icc)->info.fill,((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 137 if (((PC_Factor*)icc)->info.shiftpd) {ierr = PetscViewerASCIIPrintf(viewer," ICC: using Manteuffel shift\n");CHKERRQ(ierr);} 138 if (((PC_Factor*)icc)->info.shiftnz) {ierr = PetscViewerASCIIPrintf(viewer," ICC: using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr);} 139 if (((PC_Factor*)icc)->fact) { 140 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio needed %G\n",icc->actualfill);CHKERRQ(ierr); 141 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");CHKERRQ(ierr); 142 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 143 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 144 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 145 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 146 ierr = MatView(((PC_Factor*)icc)->fact,viewer);CHKERRQ(ierr); 147 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 148 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 149 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 150 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 151 } 152 } else if (isstring) { 153 ierr = PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)((PC_Factor*)icc)->info.levels);CHKERRQ(ierr);CHKERRQ(ierr); 154 } else { 155 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name); 156 } 157 PetscFunctionReturn(0); 158 } 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 . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default 170 - -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value 171 is optional with PETSC_TRUE being the default 172 173 Level: beginner 174 175 Concepts: incomplete Cholesky factorization 176 177 Notes: Only implemented for some matrix formats. Not implemented in parallel (for parallel use you 178 must use MATMPIROWBS, see MatCreateMPIRowbs(), this supports only ICC(0) and this is not recommended 179 unless you really want a parallel ICC). 180 181 For BAIJ matrices this implements a point block ICC. 182 183 The Manteuffel shift is only implemented for matrices with block size 1 184 185 By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftPd(pc,PETSC_FALSE); 186 to turn off the shift. 187 188 References: 189 Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST 190 http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical 191 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 192 Science and Engineering, Kluwer, pp. 167--202. 193 194 195 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 196 PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), 197 PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 198 PCFactorSetLevels(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd(), 199 200 M*/ 201 202 EXTERN_C_BEGIN 203 #undef __FUNCT__ 204 #define __FUNCT__ "PCCreate_ICC" 205 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ICC(PC pc) 206 { 207 PetscErrorCode ierr; 208 PC_ICC *icc; 209 210 PetscFunctionBegin; 211 ierr = PetscNewLog(pc,PC_ICC,&icc);CHKERRQ(ierr); 212 213 ((PC_Factor*)icc)->fact = 0; 214 ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)icc)->ordering);CHKERRQ(ierr); 215 ierr = MatFactorInfoInitialize(&((PC_Factor*)icc)->info);CHKERRQ(ierr); 216 ((PC_Factor*)icc)->info.levels = 0; 217 ((PC_Factor*)icc)->info.fill = 1.0; 218 icc->implctx = 0; 219 220 ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT; 221 ((PC_Factor*)icc)->info.shiftnz = 0.0; 222 ((PC_Factor*)icc)->info.shiftpd = 1.0; /* true */ 223 ((PC_Factor*)icc)->info.zeropivot = 1.e-12; 224 pc->data = (void*)icc; 225 226 pc->ops->apply = PCApply_ICC; 227 pc->ops->setup = PCSetup_ICC; 228 pc->ops->destroy = PCDestroy_ICC; 229 pc->ops->setfromoptions = PCSetFromOptions_ICC; 230 pc->ops->view = PCView_ICC; 231 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 232 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 233 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 234 235 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 236 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 237 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 238 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 239 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Factor", 240 PCFactorSetShiftNonzero_Factor);CHKERRQ(ierr); 241 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Factor", 242 PCFactorSetShiftPd_Factor);CHKERRQ(ierr); 243 244 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor", 245 PCFactorSetLevels_Factor);CHKERRQ(ierr); 246 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 247 PCFactorSetFill_Factor);CHKERRQ(ierr); 248 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 249 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 EXTERN_C_END 253 254 255