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