1 #define PETSCKSP_DLL 2 3 /* 4 Defines a direct factorization preconditioner for any Mat implementation 5 Note: this need not be consided a preconditioner since it supplies 6 a direct solver. 7 */ 8 #include "../src/ksp/pc/impls/factor/factor.h" /*I "petscpc.h" I*/ 9 10 typedef struct { 11 PC_Factor hdr; 12 PetscReal actualfill; /* actual fill in factor */ 13 PetscTruth inplace; /* flag indicating in-place factorization */ 14 IS row,col; /* index sets used for reordering */ 15 PetscTruth reuseordering; /* reuses previous reordering computed */ 16 PetscTruth reusefill; /* reuse fill from previous Cholesky */ 17 } PC_Cholesky; 18 19 EXTERN_C_BEGIN 20 #undef __FUNCT__ 21 #define __FUNCT__ "PCFactorSetReuseOrdering_Cholesky" 22 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering_Cholesky(PC pc,PetscTruth flag) 23 { 24 PC_Cholesky *lu; 25 26 PetscFunctionBegin; 27 lu = (PC_Cholesky*)pc->data; 28 lu->reuseordering = flag; 29 PetscFunctionReturn(0); 30 } 31 EXTERN_C_END 32 33 EXTERN_C_BEGIN 34 #undef __FUNCT__ 35 #define __FUNCT__ "PCFactorSetReuseFill_Cholesky" 36 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseFill_Cholesky(PC pc,PetscTruth flag) 37 { 38 PC_Cholesky *lu; 39 40 PetscFunctionBegin; 41 lu = (PC_Cholesky*)pc->data; 42 lu->reusefill = flag; 43 PetscFunctionReturn(0); 44 } 45 EXTERN_C_END 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "PCSetFromOptions_Cholesky" 49 static PetscErrorCode PCSetFromOptions_Cholesky(PC pc) 50 { 51 PetscErrorCode ierr; 52 53 PetscFunctionBegin; 54 ierr = PetscOptionsHead("Cholesky options");CHKERRQ(ierr); 55 ierr = PCSetFromOptions_Factor(pc);CHKERRQ(ierr); 56 ierr = PetscOptionsTail();CHKERRQ(ierr); 57 PetscFunctionReturn(0); 58 } 59 60 #undef __FUNCT__ 61 #define __FUNCT__ "PCView_Cholesky" 62 static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer) 63 { 64 PC_Cholesky *lu = (PC_Cholesky*)pc->data; 65 PetscErrorCode ierr; 66 PetscTruth iascii,isstring; 67 68 PetscFunctionBegin; 69 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 70 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 71 if (iascii) { 72 73 if (lu->inplace) {ierr = PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorization\n");CHKERRQ(ierr);} 74 else {ierr = PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorization\n");CHKERRQ(ierr);} 75 ierr = PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",((PC_Factor*)lu)->ordering);CHKERRQ(ierr); 76 if (((PC_Factor*)lu)->fact) { 77 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: factor fill ratio needed %G\n",lu->actualfill);CHKERRQ(ierr); 78 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");CHKERRQ(ierr); 79 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 80 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 81 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 82 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 83 ierr = MatView(((PC_Factor*)lu)->fact,viewer);CHKERRQ(ierr); 84 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 85 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 86 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 87 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 88 } 89 if (lu->reusefill) {ierr = PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");CHKERRQ(ierr);} 90 if (lu->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");CHKERRQ(ierr);} 91 } else if (isstring) { 92 ierr = PetscViewerStringSPrintf(viewer," order=%s",((PC_Factor*)lu)->ordering);CHKERRQ(ierr);CHKERRQ(ierr); 93 } else { 94 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCCHOLESKY",((PetscObject)viewer)->type_name); 95 } 96 PetscFunctionReturn(0); 97 } 98 99 100 #undef __FUNCT__ 101 #define __FUNCT__ "PCSetUp_Cholesky" 102 static PetscErrorCode PCSetUp_Cholesky(PC pc) 103 { 104 PetscErrorCode ierr; 105 PetscTruth flg; 106 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 107 108 PetscFunctionBegin; 109 if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill; 110 111 if (dir->inplace) { 112 if (dir->row && dir->col && (dir->row != dir->col)) { 113 ierr = ISDestroy(dir->row);CHKERRQ(ierr); 114 dir->row = 0; 115 } 116 if (dir->col) { 117 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 118 dir->col = 0; 119 } 120 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 121 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 122 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 123 dir->col=0; 124 } 125 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 126 ierr = MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 127 ((PC_Factor*)dir)->fact = pc->pmat; 128 } else { 129 MatInfo info; 130 if (!pc->setupcalled) { 131 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 132 /* check if dir->row == dir->col */ 133 ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr); 134 if (!flg) SETERRQ(PETSC_ERR_ARG_INCOMP,"row and column permutations must equal"); 135 ierr = ISDestroy(dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */ 136 dir->col=0; 137 138 flg = PETSC_FALSE; 139 ierr = PetscOptionsGetTruth(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 140 if (flg) { 141 PetscReal tol = 1.e-10; 142 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 143 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 144 } 145 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 146 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 147 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 148 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 149 dir->actualfill = info.fill_ratio_needed; 150 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 151 } else if (pc->flag != SAME_NONZERO_PATTERN) { 152 if (!dir->reuseordering) { 153 if (dir->row && dir->col && (dir->row != dir->col)) { 154 ierr = ISDestroy(dir->row);CHKERRQ(ierr); 155 dir->row = 0; 156 } 157 if (dir->col) { 158 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 159 dir->col =0; 160 } 161 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 162 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 163 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 164 dir->col=0; 165 } 166 flg = PETSC_FALSE; 167 ierr = PetscOptionsGetTruth(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 168 if (flg) { 169 PetscReal tol = 1.e-10; 170 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 171 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 172 } 173 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 174 } 175 ierr = MatDestroy(((PC_Factor*)dir)->fact);CHKERRQ(ierr); 176 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 177 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 178 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 179 dir->actualfill = info.fill_ratio_needed; 180 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 181 } 182 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 183 } 184 PetscFunctionReturn(0); 185 } 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "PCDestroy_Cholesky" 189 static PetscErrorCode PCDestroy_Cholesky(PC pc) 190 { 191 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 if (!dir->inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 196 if (dir->row) {ierr = ISDestroy(dir->row);CHKERRQ(ierr);} 197 if (dir->col) {ierr = ISDestroy(dir->col);CHKERRQ(ierr);} 198 ierr = PetscStrfree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 199 ierr = PetscStrfree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 200 ierr = PetscFree(dir);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "PCApply_Cholesky" 206 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y) 207 { 208 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 if (dir->inplace) {ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr);} 213 else {ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 214 PetscFunctionReturn(0); 215 } 216 217 #undef __FUNCT__ 218 #define __FUNCT__ "PCApplyTranspose_Cholesky" 219 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y) 220 { 221 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 222 PetscErrorCode ierr; 223 224 PetscFunctionBegin; 225 if (dir->inplace) {ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr);} 226 else {ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 227 PetscFunctionReturn(0); 228 } 229 230 /* -----------------------------------------------------------------------------------*/ 231 232 EXTERN_C_BEGIN 233 #undef __FUNCT__ 234 #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky" 235 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseInPlace_Cholesky(PC pc) 236 { 237 PC_Cholesky *dir; 238 239 PetscFunctionBegin; 240 dir = (PC_Cholesky*)pc->data; 241 dir->inplace = PETSC_TRUE; 242 PetscFunctionReturn(0); 243 } 244 EXTERN_C_END 245 246 /* -----------------------------------------------------------------------------------*/ 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "PCFactorSetReuseOrdering" 250 /*@ 251 PCFactorSetReuseOrdering - When similar matrices are factored, this 252 causes the ordering computed in the first factor to be used for all 253 following factors. 254 255 Collective on PC 256 257 Input Parameters: 258 + pc - the preconditioner context 259 - flag - PETSC_TRUE to reuse else PETSC_FALSE 260 261 Options Database Key: 262 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 263 264 Level: intermediate 265 266 .keywords: PC, levels, reordering, factorization, incomplete, LU 267 268 .seealso: PCFactorSetReuseFill() 269 @*/ 270 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering(PC pc,PetscTruth flag) 271 { 272 PetscErrorCode ierr,(*f)(PC,PetscTruth); 273 274 PetscFunctionBegin; 275 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 276 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",(void (**)(void))&f);CHKERRQ(ierr); 277 if (f) { 278 ierr = (*f)(pc,flag);CHKERRQ(ierr); 279 } 280 PetscFunctionReturn(0); 281 } 282 283 284 /*MC 285 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 286 287 Options Database Keys: 288 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 289 . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like spooles 290 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 291 . -pc_factor_fill <fill> - Sets fill amount 292 . -pc_factor_in_place - Activates in-place factorization 293 . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 294 . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default 295 - -pc_factor_shift_positive_definite true or false - Activate/Deactivate PCFactorSetShiftPd(); the value 296 is optional with true being the default 297 298 Notes: Not all options work for all matrix formats 299 300 Level: beginner 301 302 Concepts: Cholesky factorization, direct solver 303 304 Notes: Usually this will compute an "exact" solution in one iteration and does 305 not need a Krylov method (i.e. you can use -ksp_type preonly, or 306 KSPSetType(ksp,KSPPREONLY) for the Krylov method 307 308 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 309 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 310 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), PCFactorSetShiftInBlocks() 311 PCFactorSetUseInPlace(), PCFactorSetMatOrderingType() 312 313 M*/ 314 315 EXTERN_C_BEGIN 316 #undef __FUNCT__ 317 #define __FUNCT__ "PCCreate_Cholesky" 318 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Cholesky(PC pc) 319 { 320 PetscErrorCode ierr; 321 PC_Cholesky *dir; 322 323 PetscFunctionBegin; 324 ierr = PetscNewLog(pc,PC_Cholesky,&dir);CHKERRQ(ierr); 325 326 ((PC_Factor*)dir)->fact = 0; 327 dir->inplace = PETSC_FALSE; 328 ierr = MatFactorInfoInitialize(&((PC_Factor*)dir)->info);CHKERRQ(ierr); 329 ((PC_Factor*)dir)->info.fill = 5.0; 330 ((PC_Factor*)dir)->info.shiftnz = 0.0; 331 ((PC_Factor*)dir)->info.shiftpd = 0.0; /* false */ 332 ((PC_Factor*)dir)->info.pivotinblocks = 1.0; 333 dir->col = 0; 334 dir->row = 0; 335 ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 336 ierr = PetscStrallocpy(MAT_SOLVER_PETSC,&((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 337 dir->reusefill = PETSC_FALSE; 338 dir->reuseordering = PETSC_FALSE; 339 pc->data = (void*)dir; 340 341 pc->ops->destroy = PCDestroy_Cholesky; 342 pc->ops->apply = PCApply_Cholesky; 343 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 344 pc->ops->setup = PCSetUp_Cholesky; 345 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 346 pc->ops->view = PCView_Cholesky; 347 pc->ops->applyrichardson = 0; 348 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 349 350 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor", 351 PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 352 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 353 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 354 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 355 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 356 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Factor", 357 PCFactorSetShiftNonzero_Factor);CHKERRQ(ierr); 358 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Factor", 359 PCFactorSetShiftPd_Factor);CHKERRQ(ierr); 360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftInBlocks_C","PCFactorSetShiftInBlocks_Factor", 361 PCFactorSetShiftInBlocks_Factor);CHKERRQ(ierr); 362 363 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 364 PCFactorSetFill_Factor);CHKERRQ(ierr); 365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky", 366 PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr); 367 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 368 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 369 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky", 370 PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr); 371 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky", 372 PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr); 373 PetscFunctionReturn(0); 374 } 375 EXTERN_C_END 376