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