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