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