1 #define PETSCKSP_DLL 2 3 #include "src/ksp/pc/impls/factor/icc/icc.h" /*I "petscpc.h" I*/ 4 5 EXTERN_C_BEGIN 6 #undef __FUNCT__ 7 #define __FUNCT__ "PCFactorSetZeroPivot_ICC" 8 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetZeroPivot_ICC(PC pc,PetscReal z) 9 { 10 PC_ICC *icc; 11 12 PetscFunctionBegin; 13 icc = (PC_ICC*)pc->data; 14 icc->info.zeropivot = z; 15 PetscFunctionReturn(0); 16 } 17 EXTERN_C_END 18 19 EXTERN_C_BEGIN 20 #undef __FUNCT__ 21 #define __FUNCT__ "PCFactorSetShiftNonzero_ICC" 22 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftNonzero_ICC(PC pc,PetscReal shift) 23 { 24 PC_ICC *dir; 25 26 PetscFunctionBegin; 27 dir = (PC_ICC*)pc->data; 28 if (shift == (PetscReal) PETSC_DECIDE) { 29 dir->info.shiftnz = 1.e-12; 30 } else { 31 dir->info.shiftnz = shift; 32 } 33 PetscFunctionReturn(0); 34 } 35 EXTERN_C_END 36 37 EXTERN_C_BEGIN 38 #undef __FUNCT__ 39 #define __FUNCT__ "PCFactorSetShiftPd_ICC" 40 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftPd_ICC(PC pc,PetscTruth shift) 41 { 42 PC_ICC *dir; 43 44 PetscFunctionBegin; 45 dir = (PC_ICC*)pc->data; 46 if (shift) { 47 dir->info.shiftpd = 1.0; 48 } else { 49 dir->info.shiftpd = 0.0; 50 } 51 PetscFunctionReturn(0); 52 } 53 EXTERN_C_END 54 55 EXTERN_C_BEGIN 56 #undef __FUNCT__ 57 #define __FUNCT__ "PCFactorSetMatOrderingType_ICC" 58 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatOrderingType_ICC(PC pc,const MatOrderingType ordering) 59 { 60 PC_ICC *dir = (PC_ICC*)pc->data; 61 PetscErrorCode ierr; 62 63 PetscFunctionBegin; 64 ierr = PetscStrfree(dir->ordering);CHKERRQ(ierr); 65 ierr = PetscStrallocpy(ordering,&dir->ordering);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 EXTERN_C_END 69 70 EXTERN_C_BEGIN 71 #undef __FUNCT__ 72 #define __FUNCT__ "PCFactorSetFill_ICC" 73 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetFill_ICC(PC pc,PetscReal fill) 74 { 75 PC_ICC *dir; 76 77 PetscFunctionBegin; 78 dir = (PC_ICC*)pc->data; 79 dir->info.fill = fill; 80 PetscFunctionReturn(0); 81 } 82 EXTERN_C_END 83 84 EXTERN_C_BEGIN 85 #undef __FUNCT__ 86 #define __FUNCT__ "PCFactorSetLevels_ICC" 87 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetLevels_ICC(PC pc,PetscInt levels) 88 { 89 PC_ICC *icc; 90 91 PetscFunctionBegin; 92 icc = (PC_ICC*)pc->data; 93 icc->info.levels = levels; 94 PetscFunctionReturn(0); 95 } 96 EXTERN_C_END 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "PCSetup_ICC" 100 static PetscErrorCode PCSetup_ICC(PC pc) 101 { 102 PC_ICC *icc = (PC_ICC*)pc->data; 103 IS perm,cperm; 104 PetscErrorCode ierr; 105 MatInfo info; 106 107 PetscFunctionBegin; 108 ierr = MatGetOrdering(pc->pmat,icc->ordering,&perm,&cperm);CHKERRQ(ierr); 109 110 if (!pc->setupcalled) { 111 ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&icc->fact);CHKERRQ(ierr); 112 ierr = MatICCFactorSymbolic(icc->fact,pc->pmat,perm,&icc->info);CHKERRQ(ierr); 113 } else if (pc->flag != SAME_NONZERO_PATTERN) { 114 ierr = MatDestroy(icc->fact);CHKERRQ(ierr); 115 ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&icc->fact);CHKERRQ(ierr); 116 ierr = MatICCFactorSymbolic(icc->fact,pc->pmat,perm,&icc->info);CHKERRQ(ierr); 117 } 118 ierr = MatGetInfo(icc->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 119 icc->actualfill = info.fill_ratio_needed; 120 121 ierr = ISDestroy(cperm);CHKERRQ(ierr); 122 ierr = ISDestroy(perm);CHKERRQ(ierr); 123 ierr = MatCholeskyFactorNumeric(icc->fact,pc->pmat,&icc->info);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "PCDestroy_ICC" 129 static PetscErrorCode PCDestroy_ICC(PC pc) 130 { 131 PC_ICC *icc = (PC_ICC*)pc->data; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 if (icc->fact) {ierr = MatDestroy(icc->fact);CHKERRQ(ierr);} 136 ierr = PetscStrfree(icc->ordering);CHKERRQ(ierr); 137 ierr = PetscFree(icc);CHKERRQ(ierr); 138 PetscFunctionReturn(0); 139 } 140 141 #undef __FUNCT__ 142 #define __FUNCT__ "PCApply_ICC" 143 static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y) 144 { 145 PC_ICC *icc = (PC_ICC*)pc->data; 146 PetscErrorCode ierr; 147 148 PetscFunctionBegin; 149 ierr = MatSolve(icc->fact,x,y);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "PCApplySymmetricLeft_ICC" 155 static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y) 156 { 157 PetscErrorCode ierr; 158 PC_ICC *icc = (PC_ICC*)pc->data; 159 160 PetscFunctionBegin; 161 ierr = MatForwardSolve(icc->fact,x,y);CHKERRQ(ierr); 162 PetscFunctionReturn(0); 163 } 164 165 #undef __FUNCT__ 166 #define __FUNCT__ "PCApplySymmetricRight_ICC" 167 static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y) 168 { 169 PetscErrorCode ierr; 170 PC_ICC *icc = (PC_ICC*)pc->data; 171 172 PetscFunctionBegin; 173 ierr = MatBackwardSolve(icc->fact,x,y);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176 177 #undef __FUNCT__ 178 #define __FUNCT__ "PCFactorGetMatrix_ICC" 179 static PetscErrorCode PCFactorGetMatrix_ICC(PC pc,Mat *mat) 180 { 181 PC_ICC *icc = (PC_ICC*)pc->data; 182 183 PetscFunctionBegin; 184 *mat = icc->fact; 185 PetscFunctionReturn(0); 186 } 187 188 #undef __FUNCT__ 189 #define __FUNCT__ "PCSetFromOptions_ICC" 190 static PetscErrorCode PCSetFromOptions_ICC(PC pc) 191 { 192 PC_ICC *icc = (PC_ICC*)pc->data; 193 char tname[256]; 194 PetscTruth flg; 195 PetscErrorCode ierr; 196 PetscFList ordlist; 197 198 PetscFunctionBegin; 199 ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr); 200 ierr = PetscOptionsHead("ICC Options");CHKERRQ(ierr); 201 ierr = PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",icc->info.levels,&icc->info.levels,&flg);CHKERRQ(ierr); 202 ierr = PetscOptionsReal("-pc_factor_fill","Expected fill in factorization","PCFactorSetFill",icc->info.fill,&icc->info.fill,&flg);CHKERRQ(ierr); 203 ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr); 204 ierr = PetscOptionsList("-pc_factor_mat_ordering_type","Reorder to reduce nonzeros in ICC","PCFactorSetMatOrderingType",ordlist,icc->ordering,tname,256,&flg);CHKERRQ(ierr); 205 if (flg) { 206 ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr); 207 } 208 ierr = PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);CHKERRQ(ierr); 209 if (flg) { 210 ierr = PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);CHKERRQ(ierr); 211 } 212 ierr = PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",icc->info.shiftnz,&icc->info.shiftnz,0);CHKERRQ(ierr); 213 flg = (icc->info.shiftpd > 0.0) ? PETSC_TRUE : PETSC_FALSE; 214 ierr = PetscOptionsTruth("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 215 ierr = PCFactorSetShiftPd(pc,flg);CHKERRQ(ierr); 216 ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",icc->info.zeropivot,&icc->info.zeropivot,0);CHKERRQ(ierr); 217 218 ierr = PetscOptionsTail();CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 } 221 222 #undef __FUNCT__ 223 #define __FUNCT__ "PCView_ICC" 224 static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer) 225 { 226 PC_ICC *icc = (PC_ICC*)pc->data; 227 PetscErrorCode ierr; 228 PetscTruth isstring,iascii; 229 230 PetscFunctionBegin; 231 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 232 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 233 if (iascii) { 234 if (icc->info.levels == 1) { 235 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D level of fill\n",(PetscInt)icc->info.levels);CHKERRQ(ierr); 236 } else { 237 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D levels of fill\n",(PetscInt)icc->info.levels);CHKERRQ(ierr); 238 } 239 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio allocated %G, ordering used %s\n",icc->info.fill,icc->ordering);CHKERRQ(ierr); 240 if (icc->info.shiftpd) {ierr = PetscViewerASCIIPrintf(viewer," ICC: using Manteuffel shift\n");CHKERRQ(ierr);} 241 if (icc->info.shiftnz) {ierr = PetscViewerASCIIPrintf(viewer," ICC: using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr);} 242 if (icc->fact) { 243 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio needed %G\n",icc->actualfill);CHKERRQ(ierr); 244 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");CHKERRQ(ierr); 245 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 247 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 248 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 249 ierr = MatView(icc->fact,viewer);CHKERRQ(ierr); 250 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 251 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 252 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 253 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 254 } 255 } else if (isstring) { 256 ierr = PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)icc->info.levels);CHKERRQ(ierr);CHKERRQ(ierr); 257 } else { 258 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name); 259 } 260 PetscFunctionReturn(0); 261 } 262 263 /*MC 264 PCICC - Incomplete Cholesky factorization preconditioners. 265 266 Options Database Keys: 267 + -pc_factor_levels <k> - number of levels of fill for ICC(k) 268 . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for 269 its factorization (overwrites original matrix) 270 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 271 . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 272 . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default 273 - -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value 274 is optional with PETSC_TRUE being the default 275 276 Level: beginner 277 278 Concepts: incomplete Cholesky factorization 279 280 Notes: Only implemented for some matrix formats. Not implemented in parallel (for parallel use you 281 must use MATMPIROWBS, see MatCreateMPIRowbs(), this supports only ICC(0) and this is not recommended 282 unless you really want a parallel ICC). 283 284 For BAIJ matrices this implements a point block ICC. 285 286 The Manteuffel shift is only implemented for matrices with block size 1 287 288 By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftPd(pc,PETSC_FALSE); 289 to turn off the shift. 290 291 References: 292 Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST 293 http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical 294 Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in 295 Science and Engineering, Kluwer, pp. 167--202. 296 297 298 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 299 PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), 300 PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 301 PCFactorSetLevels(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd(), 302 303 M*/ 304 305 EXTERN_C_BEGIN 306 #undef __FUNCT__ 307 #define __FUNCT__ "PCCreate_ICC" 308 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ICC(PC pc) 309 { 310 PetscErrorCode ierr; 311 PC_ICC *icc; 312 313 PetscFunctionBegin; 314 ierr = PetscNewLog(pc,PC_ICC,&icc);CHKERRQ(ierr); 315 316 icc->fact = 0; 317 ierr = PetscStrallocpy(MATORDERING_NATURAL,&icc->ordering);CHKERRQ(ierr); 318 ierr = MatFactorInfoInitialize(&icc->info);CHKERRQ(ierr); 319 icc->info.levels = 0; 320 icc->info.fill = 1.0; 321 icc->implctx = 0; 322 323 icc->info.dtcol = PETSC_DEFAULT; 324 icc->info.shiftnz = 0.0; 325 icc->info.shiftpd = 1.0; /* true */ 326 icc->info.zeropivot = 1.e-12; 327 pc->data = (void*)icc; 328 329 pc->ops->apply = PCApply_ICC; 330 pc->ops->setup = PCSetup_ICC; 331 pc->ops->destroy = PCDestroy_ICC; 332 pc->ops->setfromoptions = PCSetFromOptions_ICC; 333 pc->ops->view = PCView_ICC; 334 pc->ops->getfactoredmatrix = PCFactorGetMatrix_ICC; 335 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 336 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 337 338 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_ICC", 339 PCFactorSetZeroPivot_ICC);CHKERRQ(ierr); 340 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_ICC", 341 PCFactorSetShiftNonzero_ICC);CHKERRQ(ierr); 342 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_ICC", 343 PCFactorSetShiftPd_ICC);CHKERRQ(ierr); 344 345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_ICC", 346 PCFactorSetLevels_ICC);CHKERRQ(ierr); 347 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_ICC", 348 PCFactorSetFill_ICC);CHKERRQ(ierr); 349 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_ICC", 350 PCFactorSetMatOrderingType_ICC);CHKERRQ(ierr); 351 PetscFunctionReturn(0); 352 } 353 EXTERN_C_END 354 355 356