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__ "PCGetFactoredMatrix_ICC" 183 static PetscErrorCode PCGetFactoredMatrix_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 ierr = PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShift",&flg);CHKERRQ(ierr); 218 if (flg) { 219 ierr = PCFactorSetShiftPd(pc,PETSC_TRUE);CHKERRQ(ierr); 220 } else { 221 ierr = PCFactorSetShiftPd(pc,PETSC_FALSE);CHKERRQ(ierr); 222 } 223 ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",icc->info.zeropivot,&icc->info.zeropivot,0);CHKERRQ(ierr); 224 225 ierr = PetscOptionsTail();CHKERRQ(ierr); 226 PetscFunctionReturn(0); 227 } 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "PCView_ICC" 231 static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer) 232 { 233 PC_ICC *icc = (PC_ICC*)pc->data; 234 PetscErrorCode ierr; 235 PetscTruth isstring,iascii; 236 237 PetscFunctionBegin; 238 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 239 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 240 if (iascii) { 241 if (icc->info.levels == 1) { 242 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D level of fill\n",(PetscInt)icc->info.levels);CHKERRQ(ierr); 243 } else { 244 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D levels of fill\n",(PetscInt)icc->info.levels);CHKERRQ(ierr); 245 } 246 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio allocated %G\n",icc->info.fill);CHKERRQ(ierr); 247 if (icc->info.shiftpd) {ierr = PetscViewerASCIIPrintf(viewer," ICC: using Manteuffel shift\n");CHKERRQ(ierr);} 248 if (icc->fact) { 249 ierr = PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio needed %G\n",icc->actualfill);CHKERRQ(ierr); 250 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");CHKERRQ(ierr); 251 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 252 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 253 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 254 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 255 ierr = MatView(icc->fact,viewer);CHKERRQ(ierr); 256 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 257 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 259 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 260 } 261 } else if (isstring) { 262 ierr = PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)icc->info.levels);CHKERRQ(ierr);CHKERRQ(ierr); 263 } else { 264 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name); 265 } 266 PetscFunctionReturn(0); 267 } 268 269 /*MC 270 PCICC - Incomplete Cholesky factorization preconditioners. 271 272 Options Database Keys: 273 + -pc_factor_levels <k> - number of levels of fill for ICC(k) 274 . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for 275 its factorization (overwrites original matrix) 276 . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 277 . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 278 . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default 279 - -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value 280 is optional with PETSC_TRUE being the default 281 282 Level: beginner 283 284 Concepts: incomplete Cholesky factorization 285 286 Notes: Only implemented for some matrix formats. Not implemented in parallel 287 288 For BAIJ matrices this implements a point block ICC. 289 290 The Manteuffel shift is only implemented for matrices with block size 1 291 292 By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftPd(pc,PETSC_FALSE); 293 to turn off the shift. 294 295 296 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 297 PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), 298 PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(), 299 PCFactorSetLevels(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd(), 300 301 M*/ 302 303 EXTERN_C_BEGIN 304 #undef __FUNCT__ 305 #define __FUNCT__ "PCCreate_ICC" 306 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ICC(PC pc) 307 { 308 PetscErrorCode ierr; 309 PC_ICC *icc; 310 311 PetscFunctionBegin; 312 ierr = PetscNewLog(pc,PC_ICC,&icc);CHKERRQ(ierr); 313 314 icc->fact = 0; 315 ierr = PetscStrallocpy(MATORDERING_NATURAL,&icc->ordering);CHKERRQ(ierr); 316 ierr = MatFactorInfoInitialize(&icc->info);CHKERRQ(ierr); 317 icc->info.levels = 0; 318 icc->info.fill = 1.0; 319 icc->implctx = 0; 320 321 icc->info.dtcol = PETSC_DEFAULT; 322 icc->info.shiftnz = 0.0; 323 icc->info.shiftpd = 1.0; /* true */ 324 icc->info.shift_fraction = 0.0; 325 icc->info.zeropivot = 1.e-12; 326 pc->data = (void*)icc; 327 328 pc->ops->apply = PCApply_ICC; 329 pc->ops->setup = PCSetup_ICC; 330 pc->ops->destroy = PCDestroy_ICC; 331 pc->ops->setfromoptions = PCSetFromOptions_ICC; 332 pc->ops->view = PCView_ICC; 333 pc->ops->getfactoredmatrix = PCGetFactoredMatrix_ICC; 334 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 335 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 336 337 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_ICC", 338 PCFactorSetZeroPivot_ICC);CHKERRQ(ierr); 339 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_ICC", 340 PCFactorSetShiftNonzero_ICC);CHKERRQ(ierr); 341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_ICC", 342 PCFactorSetShiftPd_ICC);CHKERRQ(ierr); 343 344 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_ICC", 345 PCFactorSetLevels_ICC);CHKERRQ(ierr); 346 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_ICC", 347 PCFactorSetFill_ICC);CHKERRQ(ierr); 348 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_ICC", 349 PCFactorSetMatOrderingType_ICC);CHKERRQ(ierr); 350 PetscFunctionReturn(0); 351 } 352 EXTERN_C_END 353 354 355