1 /* 2 Defines a Cholesky factorization preconditioner for any Mat implementation. 3 Presently only provided for MPIRowbs format (i.e. BlockSolve). 4 */ 5 6 #include "src/ksp/pc/impls/factor/icc/icc.h" /*I "petscpc.h" I*/ 7 8 EXTERN_C_BEGIN 9 #undef __FUNCT__ 10 #define __FUNCT__ "PCFactorSetZeroPivot_ICC" 11 PetscErrorCode PCFactorSetZeroPivot_ICC(PC pc,PetscReal z) 12 { 13 PC_ICC *icc; 14 15 PetscFunctionBegin; 16 icc = (PC_ICC*)pc->data; 17 icc->info.zeropivot = z; 18 PetscFunctionReturn(0); 19 } 20 EXTERN_C_END 21 22 EXTERN_C_BEGIN 23 #undef __FUNCT__ 24 #define __FUNCT__ "PCFactorSetShiftNonzero_ICC" 25 PetscErrorCode PCFactorSetShiftNonzero_ICC(PC pc,PetscReal shift) 26 { 27 PC_ICC *dir; 28 29 PetscFunctionBegin; 30 dir = (PC_ICC*)pc->data; 31 if (shift == (PetscReal) PETSC_DECIDE) { 32 dir->info.shiftnz = 1.e-12; 33 } else { 34 dir->info.shiftnz = shift; 35 } 36 PetscFunctionReturn(0); 37 } 38 EXTERN_C_END 39 40 EXTERN_C_BEGIN 41 #undef __FUNCT__ 42 #define __FUNCT__ "PCFactorSetShiftPd_ICC" 43 PetscErrorCode PCFactorSetShiftPd_ICC(PC pc,PetscTruth shift) 44 { 45 PC_ICC *dir; 46 47 PetscFunctionBegin; 48 dir = (PC_ICC*)pc->data; 49 dir->info.shiftpd = shift; 50 if (shift) dir->info.shift_fraction = 0.0; 51 PetscFunctionReturn(0); 52 } 53 EXTERN_C_END 54 55 EXTERN_C_BEGIN 56 #undef __FUNCT__ 57 #define __FUNCT__ "PCICCSetMatOrdering_ICC" 58 PetscErrorCode PCICCSetMatOrdering_ICC(PC pc,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__ "PCICCSetFill_ICC" 73 PetscErrorCode PCICCSetFill_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__ "PCICCSetLevels_ICC" 87 PetscErrorCode PCICCSetLevels_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__ "PCICCSetMatOrdering" 100 /*@ 101 PCICCSetMatOrdering - Sets the ordering routine (to reduce fill) to 102 be used it the ICC factorization. 103 104 Collective on PC 105 106 Input Parameters: 107 + pc - the preconditioner context 108 - ordering - the matrix ordering name, for example, MATORDERING_ND or MATORDERING_RCM 109 110 Options Database Key: 111 . -pc_icc_mat_ordering_type <nd,rcm,...> - Sets ordering routine 112 113 Level: intermediate 114 115 .seealso: PCLUSetMatOrdering() 116 117 .keywords: PC, ICC, set, matrix, reordering 118 119 @*/ 120 PetscErrorCode PCICCSetMatOrdering(PC pc,MatOrderingType ordering) 121 { 122 PetscErrorCode ierr,(*f)(PC,MatOrderingType); 123 124 PetscFunctionBegin; 125 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 126 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCICCSetMatOrdering_C",(void (**)(void))&f);CHKERRQ(ierr); 127 if (f) { 128 ierr = (*f)(pc,ordering);CHKERRQ(ierr); 129 } 130 PetscFunctionReturn(0); 131 } 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "PCICCSetLevels" 135 /*@ 136 PCICCSetLevels - Sets the number of levels of fill to use. 137 138 Collective on PC 139 140 Input Parameters: 141 + pc - the preconditioner context 142 - levels - number of levels of fill 143 144 Options Database Key: 145 . -pc_icc_levels <levels> - Sets fill level 146 147 Level: intermediate 148 149 Concepts: ICC^setting levels of fill 150 151 @*/ 152 PetscErrorCode PCICCSetLevels(PC pc,PetscInt levels) 153 { 154 PetscErrorCode ierr,(*f)(PC,PetscInt); 155 156 PetscFunctionBegin; 157 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 158 if (levels < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"negative levels"); 159 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCICCSetLevels_C",(void (**)(void))&f);CHKERRQ(ierr); 160 if (f) { 161 ierr = (*f)(pc,levels);CHKERRQ(ierr); 162 } 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "PCICCSetFill" 168 /*@ 169 PCICCSetFill - Indicate the amount of fill you expect in the factored matrix, 170 where fill = number nonzeros in factor/number nonzeros in original matrix. 171 172 Collective on PC 173 174 Input Parameters: 175 + pc - the preconditioner context 176 - fill - amount of expected fill 177 178 Options Database Key: 179 $ -pc_icc_fill <fill> 180 181 Note: 182 For sparse matrix factorizations it is difficult to predict how much 183 fill to expect. By running with the option -log_info PETSc will print the 184 actual amount of fill used; allowing you to set the value accurately for 185 future runs. But default PETSc uses a value of 1.0 186 187 Level: intermediate 188 189 .keywords: PC, set, factorization, direct, fill 190 191 .seealso: PCLUSetFill() 192 @*/ 193 PetscErrorCode PCICCSetFill(PC pc,PetscReal fill) 194 { 195 PetscErrorCode ierr,(*f)(PC,PetscReal); 196 197 PetscFunctionBegin; 198 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 199 if (fill < 1.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Fill factor cannot be less than 1.0"); 200 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCICCSetFill_C",(void (**)(void))&f);CHKERRQ(ierr); 201 if (f) { 202 ierr = (*f)(pc,fill);CHKERRQ(ierr); 203 } 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "PCSetup_ICC" 209 static PetscErrorCode PCSetup_ICC(PC pc) 210 { 211 PC_ICC *icc = (PC_ICC*)pc->data; 212 IS perm,cperm; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 ierr = MatGetOrdering(pc->pmat,icc->ordering,&perm,&cperm);CHKERRQ(ierr); 217 218 if (!pc->setupcalled) { 219 ierr = MatICCFactorSymbolic(pc->pmat,perm,&icc->info,&icc->fact);CHKERRQ(ierr); 220 } else if (pc->flag != SAME_NONZERO_PATTERN) { 221 ierr = MatDestroy(icc->fact);CHKERRQ(ierr); 222 ierr = MatICCFactorSymbolic(pc->pmat,perm,&icc->info,&icc->fact);CHKERRQ(ierr); 223 } 224 ierr = ISDestroy(cperm);CHKERRQ(ierr); 225 ierr = ISDestroy(perm);CHKERRQ(ierr); 226 ierr = MatCholeskyFactorNumeric(pc->pmat,&icc->info,&icc->fact);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 #undef __FUNCT__ 231 #define __FUNCT__ "PCDestroy_ICC" 232 static PetscErrorCode PCDestroy_ICC(PC pc) 233 { 234 PC_ICC *icc = (PC_ICC*)pc->data; 235 PetscErrorCode ierr; 236 237 PetscFunctionBegin; 238 if (icc->fact) {ierr = MatDestroy(icc->fact);CHKERRQ(ierr);} 239 ierr = PetscStrfree(icc->ordering);CHKERRQ(ierr); 240 ierr = PetscFree(icc);CHKERRQ(ierr); 241 PetscFunctionReturn(0); 242 } 243 244 #undef __FUNCT__ 245 #define __FUNCT__ "PCApply_ICC" 246 static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y) 247 { 248 PC_ICC *icc = (PC_ICC*)pc->data; 249 PetscErrorCode ierr; 250 251 PetscFunctionBegin; 252 ierr = MatSolve(icc->fact,x,y);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 #undef __FUNCT__ 257 #define __FUNCT__ "PCApplySymmetricLeft_ICC" 258 static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y) 259 { 260 PetscErrorCode ierr; 261 PC_ICC *icc = (PC_ICC*)pc->data; 262 263 PetscFunctionBegin; 264 ierr = MatForwardSolve(icc->fact,x,y);CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "PCApplySymmetricRight_ICC" 270 static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y) 271 { 272 PetscErrorCode ierr; 273 PC_ICC *icc = (PC_ICC*)pc->data; 274 275 PetscFunctionBegin; 276 ierr = MatBackwardSolve(icc->fact,x,y);CHKERRQ(ierr); 277 PetscFunctionReturn(0); 278 } 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "PCGetFactoredMatrix_ICC" 282 static PetscErrorCode PCGetFactoredMatrix_ICC(PC pc,Mat *mat) 283 { 284 PC_ICC *icc = (PC_ICC*)pc->data; 285 286 PetscFunctionBegin; 287 *mat = icc->fact; 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "PCSetFromOptions_ICC" 293 static PetscErrorCode PCSetFromOptions_ICC(PC pc) 294 { 295 PC_ICC *icc = (PC_ICC*)pc->data; 296 char tname[256]; 297 PetscTruth flg; 298 PetscErrorCode ierr; 299 PetscFList ordlist; 300 301 PetscFunctionBegin; 302 ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr); 303 ierr = PetscOptionsHead("ICC Options");CHKERRQ(ierr); 304 ierr = PetscOptionsReal("-pc_icc_levels","levels of fill","PCICCSetLevels",icc->info.levels,&icc->info.levels,&flg);CHKERRQ(ierr); 305 ierr = PetscOptionsReal("-pc_icc_fill","Expected fill in factorization","PCICCSetFill",icc->info.fill,&icc->info.fill,&flg);CHKERRQ(ierr); 306 ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr); 307 ierr = PetscOptionsList("-pc_icc_mat_ordering_type","Reorder to reduce nonzeros in ICC","PCICCSetMatOrdering",ordlist,icc->ordering,tname,256,&flg);CHKERRQ(ierr); 308 if (flg) { 309 ierr = PCICCSetMatOrdering(pc,tname);CHKERRQ(ierr); 310 } 311 ierr = PetscOptionsName("-pc_factor_shiftnonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);CHKERRQ(ierr); 312 if (flg) { 313 ierr = PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);CHKERRQ(ierr); 314 } 315 ierr = PetscOptionsReal("-pc_factor_shiftnonzero","Shift added to diagonal","PCFactorSetShiftNonzero",icc->info.shiftnz,&icc->info.shiftnz,0);CHKERRQ(ierr); 316 ierr = PetscOptionsName("-pc_factor_shiftpd","Manteuffel shift applied to diagonal","PCICCSetShift",&flg);CHKERRQ(ierr); 317 if (flg) { 318 ierr = PCFactorSetShiftPd(pc,PETSC_TRUE);CHKERRQ(ierr); 319 } else { 320 ierr = PCFactorSetShiftPd(pc,PETSC_FALSE);CHKERRQ(ierr); 321 } 322 ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",icc->info.zeropivot,&icc->info.zeropivot,0);CHKERRQ(ierr); 323 324 ierr = PetscOptionsTail();CHKERRQ(ierr); 325 PetscFunctionReturn(0); 326 } 327 328 #undef __FUNCT__ 329 #define __FUNCT__ "PCView_ICC" 330 static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer) 331 { 332 PC_ICC *icc = (PC_ICC*)pc->data; 333 PetscErrorCode ierr; 334 PetscTruth isstring,iascii; 335 336 PetscFunctionBegin; 337 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 338 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 339 if (iascii) { 340 if (icc->info.levels == 1) { 341 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D level of fill\n",(PetscInt)icc->info.levels);CHKERRQ(ierr); 342 } else { 343 ierr = PetscViewerASCIIPrintf(viewer," ICC: %D levels of fill\n",(PetscInt)icc->info.levels);CHKERRQ(ierr); 344 } 345 ierr = PetscViewerASCIIPrintf(viewer," ICC: max fill ratio allocated %g\n",icc->info.fill);CHKERRQ(ierr); 346 if (icc->info.shiftpd) {ierr = PetscViewerASCIIPrintf(viewer," ICC: using Manteuffel shift\n");CHKERRQ(ierr);} 347 } else if (isstring) { 348 ierr = PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)icc->info.levels);CHKERRQ(ierr);CHKERRQ(ierr); 349 } else { 350 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name); 351 } 352 PetscFunctionReturn(0); 353 } 354 355 /*MC 356 PCICC - Incomplete Cholesky factorization preconditioners. 357 358 Options Database Keys: 359 + -pc_icc_levels <k> - number of levels of fill for ICC(k) 360 . -pc_icc_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for 361 its factorization (overwrites original matrix) 362 . -pc_icc_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1 363 - -pc_icc_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix 364 365 Level: beginner 366 367 Concepts: incomplete Cholesky factorization 368 369 Notes: Only implemented for some matrix formats. Not implemented in parallel 370 371 For BAIJ matrices this implements a point block ICC. 372 373 The Manteuffel shift is only implemented for matrices with block size 1 374 375 By default, the Manteuffel is applied (for matrices with block size 1). Call PCICCSetShift(pc,PETSC_FALSE); 376 to turn off the shift. 377 378 379 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType, 380 PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), 381 PCICCSetFill(), PCICCSetMatOrdering(), PCICCSetReuseOrdering(), 382 PCICCSetLevels() 383 384 M*/ 385 386 EXTERN_C_BEGIN 387 #undef __FUNCT__ 388 #define __FUNCT__ "PCCreate_ICC" 389 PetscErrorCode PCCreate_ICC(PC pc) 390 { 391 PetscErrorCode ierr; 392 PC_ICC *icc; 393 394 PetscFunctionBegin; 395 ierr = PetscNew(PC_ICC,&icc);CHKERRQ(ierr); 396 PetscLogObjectMemory(pc,sizeof(PC_ICC)); 397 398 icc->fact = 0; 399 ierr = PetscStrallocpy(MATORDERING_NATURAL,&icc->ordering);CHKERRQ(ierr); 400 ierr = MatFactorInfoInitialize(&icc->info);CHKERRQ(ierr); 401 icc->info.levels = 0; 402 icc->info.fill = 1.0; 403 icc->implctx = 0; 404 405 icc->info.dtcol = PETSC_DEFAULT; 406 icc->info.shiftnz = 0.0; 407 icc->info.shiftpd = PETSC_TRUE; 408 icc->info.shift_fraction = 0.0; 409 icc->info.zeropivot = 1.e-12; 410 pc->data = (void*)icc; 411 412 pc->ops->apply = PCApply_ICC; 413 pc->ops->setup = PCSetup_ICC; 414 pc->ops->destroy = PCDestroy_ICC; 415 pc->ops->setfromoptions = PCSetFromOptions_ICC; 416 pc->ops->view = PCView_ICC; 417 pc->ops->getfactoredmatrix = PCGetFactoredMatrix_ICC; 418 pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC; 419 pc->ops->applysymmetricright = PCApplySymmetricRight_ICC; 420 421 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_ICC", 422 PCFactorSetZeroPivot_ICC);CHKERRQ(ierr); 423 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_ICC", 424 PCFactorSetShiftNonzero_ICC);CHKERRQ(ierr); 425 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_ICC", 426 PCFactorSetShiftPd_ICC);CHKERRQ(ierr); 427 428 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCICCSetLevels_C","PCICCSetLevels_ICC", 429 PCICCSetLevels_ICC);CHKERRQ(ierr); 430 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCICCSetFill_C","PCICCSetFill_ICC", 431 PCICCSetFill_ICC);CHKERRQ(ierr); 432 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCICCSetMatOrdering_C","PCICCSetMatOrdering_ICC", 433 PCICCSetMatOrdering_ICC);CHKERRQ(ierr); 434 PetscFunctionReturn(0); 435 } 436 EXTERN_C_END 437 438 439