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 *chol = (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 if (((PC_Factor*)chol)->info.shifttype==MAT_SHIFT_POSITIVE_DEFINITE) { 73 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: using Manteuffel shift\n");CHKERRQ(ierr); 74 } 75 if (((PC_Factor*)chol)->info.shifttype==MAT_SHIFT_NONZERO) { 76 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr); 77 } 78 if (((PC_Factor*)chol)->info.shifttype==MAT_SHIFT_INBLOCKS) { 79 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: using diagonal shift on blocks to prevent zero pivot\n");CHKERRQ(ierr); 80 } 81 82 if (chol->inplace) {ierr = PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorization\n");CHKERRQ(ierr);} 83 else {ierr = PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorization\n");CHKERRQ(ierr);} 84 ierr = PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",((PC_Factor*)chol)->ordering);CHKERRQ(ierr); 85 if (((PC_Factor*)chol)->fact) { 86 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: factor fill ratio needed %G\n",chol->actualfill);CHKERRQ(ierr); 87 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");CHKERRQ(ierr); 88 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 89 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 90 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 91 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 92 ierr = MatView(((PC_Factor*)chol)->fact,viewer);CHKERRQ(ierr); 93 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 94 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 95 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 96 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 97 } 98 if (chol->reusefill) {ierr = PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");CHKERRQ(ierr);} 99 if (chol->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");CHKERRQ(ierr);} 100 } else if (isstring) { 101 ierr = PetscViewerStringSPrintf(viewer," order=%s",((PC_Factor*)chol)->ordering);CHKERRQ(ierr);CHKERRQ(ierr); 102 } else { 103 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCCHOLESKY",((PetscObject)viewer)->type_name); 104 } 105 PetscFunctionReturn(0); 106 } 107 108 109 #undef __FUNCT__ 110 #define __FUNCT__ "PCSetUp_Cholesky" 111 static PetscErrorCode PCSetUp_Cholesky(PC pc) 112 { 113 PetscErrorCode ierr; 114 PetscTruth flg; 115 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 116 117 PetscFunctionBegin; 118 if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill; 119 120 if (dir->inplace) { 121 if (dir->row && dir->col && (dir->row != dir->col)) { 122 ierr = ISDestroy(dir->row);CHKERRQ(ierr); 123 dir->row = 0; 124 } 125 if (dir->col) { 126 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 127 dir->col = 0; 128 } 129 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 130 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 131 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 132 dir->col=0; 133 } 134 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 135 ierr = MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 136 ((PC_Factor*)dir)->fact = pc->pmat; 137 } else { 138 MatInfo info; 139 if (!pc->setupcalled) { 140 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 141 /* check if dir->row == dir->col */ 142 ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr); 143 if (!flg) SETERRQ(PETSC_ERR_ARG_INCOMP,"row and column permutations must equal"); 144 ierr = ISDestroy(dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */ 145 dir->col=0; 146 147 flg = PETSC_FALSE; 148 ierr = PetscOptionsGetTruth(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 149 if (flg) { 150 PetscReal tol = 1.e-10; 151 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 152 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 153 } 154 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 155 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 156 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 157 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 158 dir->actualfill = info.fill_ratio_needed; 159 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 160 } else if (pc->flag != SAME_NONZERO_PATTERN) { 161 if (!dir->reuseordering) { 162 if (dir->row && dir->col && (dir->row != dir->col)) { 163 ierr = ISDestroy(dir->row);CHKERRQ(ierr); 164 dir->row = 0; 165 } 166 if (dir->col) { 167 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 168 dir->col =0; 169 } 170 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 171 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 172 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 173 dir->col=0; 174 } 175 flg = PETSC_FALSE; 176 ierr = PetscOptionsGetTruth(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 177 if (flg) { 178 PetscReal tol = 1.e-10; 179 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 180 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 181 } 182 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 183 } 184 ierr = MatDestroy(((PC_Factor*)dir)->fact);CHKERRQ(ierr); 185 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 186 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 187 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 188 dir->actualfill = info.fill_ratio_needed; 189 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 190 } 191 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 192 } 193 PetscFunctionReturn(0); 194 } 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "PCDestroy_Cholesky" 198 static PetscErrorCode PCDestroy_Cholesky(PC pc) 199 { 200 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 201 PetscErrorCode ierr; 202 203 PetscFunctionBegin; 204 if (!dir->inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 205 if (dir->row) {ierr = ISDestroy(dir->row);CHKERRQ(ierr);} 206 if (dir->col) {ierr = ISDestroy(dir->col);CHKERRQ(ierr);} 207 ierr = PetscStrfree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 208 ierr = PetscStrfree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 209 ierr = PetscFree(dir);CHKERRQ(ierr); 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "PCApply_Cholesky" 215 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y) 216 { 217 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 218 PetscErrorCode ierr; 219 220 PetscFunctionBegin; 221 if (dir->inplace) {ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr);} 222 else {ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 223 PetscFunctionReturn(0); 224 } 225 226 #undef __FUNCT__ 227 #define __FUNCT__ "PCApplyTranspose_Cholesky" 228 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y) 229 { 230 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 231 PetscErrorCode ierr; 232 233 PetscFunctionBegin; 234 if (dir->inplace) {ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr);} 235 else {ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 236 PetscFunctionReturn(0); 237 } 238 239 /* -----------------------------------------------------------------------------------*/ 240 241 EXTERN_C_BEGIN 242 #undef __FUNCT__ 243 #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky" 244 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseInPlace_Cholesky(PC pc) 245 { 246 PC_Cholesky *dir; 247 248 PetscFunctionBegin; 249 dir = (PC_Cholesky*)pc->data; 250 dir->inplace = PETSC_TRUE; 251 PetscFunctionReturn(0); 252 } 253 EXTERN_C_END 254 255 /* -----------------------------------------------------------------------------------*/ 256 257 #undef __FUNCT__ 258 #define __FUNCT__ "PCFactorSetReuseOrdering" 259 /*@ 260 PCFactorSetReuseOrdering - When similar matrices are factored, this 261 causes the ordering computed in the first factor to be used for all 262 following factors. 263 264 Collective on PC 265 266 Input Parameters: 267 + pc - the preconditioner context 268 - flag - PETSC_TRUE to reuse else PETSC_FALSE 269 270 Options Database Key: 271 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 272 273 Level: intermediate 274 275 .keywords: PC, levels, reordering, factorization, incomplete, LU 276 277 .seealso: PCFactorSetReuseFill() 278 @*/ 279 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering(PC pc,PetscTruth flag) 280 { 281 PetscErrorCode ierr,(*f)(PC,PetscTruth); 282 283 PetscFunctionBegin; 284 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 285 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",(void (**)(void))&f);CHKERRQ(ierr); 286 if (f) { 287 ierr = (*f)(pc,flag);CHKERRQ(ierr); 288 } 289 PetscFunctionReturn(0); 290 } 291 292 293 /*MC 294 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 295 296 Options Database Keys: 297 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 298 . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like spooles 299 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 300 . -pc_factor_fill <fill> - Sets fill amount 301 . -pc_factor_in_place - Activates in-place factorization 302 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 303 304 Notes: Not all options work for all matrix formats 305 306 Level: beginner 307 308 Concepts: Cholesky factorization, direct solver 309 310 Notes: Usually this will compute an "exact" solution in one iteration and does 311 not need a Krylov method (i.e. you can use -ksp_type preonly, or 312 KSPSetType(ksp,KSPPREONLY) for the Krylov method 313 314 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 315 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 316 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 317 PCFactorSetUseInPlace(), PCFactorSetMatOrderingType() 318 319 M*/ 320 321 EXTERN_C_BEGIN 322 #undef __FUNCT__ 323 #define __FUNCT__ "PCCreate_Cholesky" 324 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Cholesky(PC pc) 325 { 326 PetscErrorCode ierr; 327 PC_Cholesky *dir; 328 329 PetscFunctionBegin; 330 ierr = PetscNewLog(pc,PC_Cholesky,&dir);CHKERRQ(ierr); 331 332 ((PC_Factor*)dir)->fact = 0; 333 dir->inplace = PETSC_FALSE; 334 ierr = MatFactorInfoInitialize(&((PC_Factor*)dir)->info);CHKERRQ(ierr); 335 ((PC_Factor*)dir)->info.fill = 5.0; 336 ((PC_Factor*)dir)->info.shifttype = MAT_SHIFT_NONE; 337 ((PC_Factor*)dir)->info.shiftamount = 0.0; 338 ((PC_Factor*)dir)->info.pivotinblocks = 1.0; 339 dir->col = 0; 340 dir->row = 0; 341 ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 342 ierr = PetscStrallocpy(MAT_SOLVER_PETSC,&((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 343 dir->reusefill = PETSC_FALSE; 344 dir->reuseordering = PETSC_FALSE; 345 pc->data = (void*)dir; 346 347 pc->ops->destroy = PCDestroy_Cholesky; 348 pc->ops->apply = PCApply_Cholesky; 349 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 350 pc->ops->setup = PCSetUp_Cholesky; 351 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 352 pc->ops->view = PCView_Cholesky; 353 pc->ops->applyrichardson = 0; 354 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 355 356 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor", 357 PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 358 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 359 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 361 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 362 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor", 363 PCFactorSetShiftType_Factor);CHKERRQ(ierr); 364 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor", 365 PCFactorSetShiftAmount_Factor);CHKERRQ(ierr); 366 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 367 PCFactorSetFill_Factor);CHKERRQ(ierr); 368 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky", 369 PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr); 370 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 371 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 372 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky", 373 PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr); 374 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky", 375 PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr); 376 PetscFunctionReturn(0); 377 } 378 EXTERN_C_END 379