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