1 2 /* 3 Defines a direct factorization preconditioner for any Mat implementation 4 Note: this need not be consided a preconditioner since it supplies 5 a direct solver. 6 */ 7 #include <../src/ksp/pc/impls/factor/factor.h> /*I "petscpc.h" I*/ 8 9 typedef struct { 10 PC_Factor hdr; 11 IS row, col; /* index sets used for reordering */ 12 } PC_Cholesky; 13 14 static PetscErrorCode PCSetFromOptions_Cholesky(PC pc, PetscOptionItems *PetscOptionsObject) { 15 PetscFunctionBegin; 16 PetscOptionsHeadBegin(PetscOptionsObject, "Cholesky options"); 17 PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject)); 18 PetscOptionsHeadEnd(); 19 PetscFunctionReturn(0); 20 } 21 22 static PetscErrorCode PCSetUp_Cholesky(PC pc) { 23 PetscBool flg; 24 PC_Cholesky *dir = (PC_Cholesky *)pc->data; 25 MatSolverType stype; 26 MatFactorError err; 27 const char *prefix; 28 29 PetscFunctionBegin; 30 pc->failedreason = PC_NOERROR; 31 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill; 32 33 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 34 PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix)); 35 36 PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure)); 37 if (dir->hdr.inplace) { 38 if (dir->row && dir->col && (dir->row != dir->col)) PetscCall(ISDestroy(&dir->row)); 39 PetscCall(ISDestroy(&dir->col)); 40 /* should only get reordering if the factor matrix uses it but cannot determine because MatGetFactor() not called */ 41 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 42 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col)); 43 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 44 PetscCall(ISDestroy(&dir->col)); 45 } 46 if (dir->row) PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)dir->row)); 47 PetscCall(MatCholeskyFactor(pc->pmat, dir->row, &((PC_Factor *)dir)->info)); 48 PetscCall(MatFactorGetError(pc->pmat, &err)); 49 if (err) { /* Factor() fails */ 50 pc->failedreason = (PCFailedReason)err; 51 PetscFunctionReturn(0); 52 } 53 54 ((PC_Factor *)dir)->fact = pc->pmat; 55 } else { 56 MatInfo info; 57 58 if (!pc->setupcalled) { 59 PetscBool canuseordering; 60 if (!((PC_Factor *)dir)->fact) { 61 PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_CHOLESKY, &((PC_Factor *)dir)->fact)); 62 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)((PC_Factor *)dir)->fact)); 63 } 64 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering)); 65 if (canuseordering) { 66 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 67 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col)); 68 /* check if dir->row == dir->col */ 69 if (dir->row) { 70 PetscCall(ISEqual(dir->row, dir->col, &flg)); 71 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "row and column permutations must be equal"); 72 } 73 PetscCall(ISDestroy(&dir->col)); /* only pass one ordering into CholeskyFactor */ 74 75 flg = PETSC_FALSE; 76 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg, NULL)); 77 if (flg) { 78 PetscReal tol = 1.e-10; 79 PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL)); 80 PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row)); 81 } 82 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)dir->row)); 83 } 84 PetscCall(MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info)); 85 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info)); 86 dir->hdr.actualfill = info.fill_ratio_needed; 87 } else if (pc->flag != SAME_NONZERO_PATTERN) { 88 if (!dir->hdr.reuseordering) { 89 PetscBool canuseordering; 90 PetscCall(MatDestroy(&((PC_Factor *)dir)->fact)); 91 PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_CHOLESKY, &((PC_Factor *)dir)->fact)); 92 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)((PC_Factor *)dir)->fact)); 93 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering)); 94 if (canuseordering) { 95 PetscCall(ISDestroy(&dir->row)); 96 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 97 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col)); 98 PetscCall(ISDestroy(&dir->col)); /* only use dir->row ordering in CholeskyFactor */ 99 100 flg = PETSC_FALSE; 101 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg, NULL)); 102 if (flg) { 103 PetscReal tol = 1.e-10; 104 PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL)); 105 PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row)); 106 } 107 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)dir->row)); 108 } 109 } 110 PetscCall(MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info)); 111 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info)); 112 dir->hdr.actualfill = info.fill_ratio_needed; 113 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)((PC_Factor *)dir)->fact)); 114 } else { 115 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 116 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) { 117 PetscCall(MatFactorClearError(((PC_Factor *)dir)->fact)); 118 pc->failedreason = PC_NOERROR; 119 } 120 } 121 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 122 if (err) { /* FactorSymbolic() fails */ 123 pc->failedreason = (PCFailedReason)err; 124 PetscFunctionReturn(0); 125 } 126 127 PetscCall(MatCholeskyFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info)); 128 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 129 if (err) { /* FactorNumeric() fails */ 130 pc->failedreason = (PCFailedReason)err; 131 } 132 } 133 134 PetscCall(PCFactorGetMatSolverType(pc, &stype)); 135 if (!stype) { 136 MatSolverType solverpackage; 137 PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage)); 138 PetscCall(PCFactorSetMatSolverType(pc, solverpackage)); 139 } 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode PCReset_Cholesky(PC pc) { 144 PC_Cholesky *dir = (PC_Cholesky *)pc->data; 145 146 PetscFunctionBegin; 147 if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact)); 148 PetscCall(ISDestroy(&dir->row)); 149 PetscCall(ISDestroy(&dir->col)); 150 PetscFunctionReturn(0); 151 } 152 153 static PetscErrorCode PCDestroy_Cholesky(PC pc) { 154 PC_Cholesky *dir = (PC_Cholesky *)pc->data; 155 156 PetscFunctionBegin; 157 PetscCall(PCReset_Cholesky(pc)); 158 PetscCall(PetscFree(((PC_Factor *)dir)->ordering)); 159 PetscCall(PetscFree(((PC_Factor *)dir)->solvertype)); 160 PetscCall(PCFactorClearComposedFunctions(pc)); 161 PetscCall(PetscFree(pc->data)); 162 PetscFunctionReturn(0); 163 } 164 165 static PetscErrorCode PCApply_Cholesky(PC pc, Vec x, Vec y) { 166 PC_Cholesky *dir = (PC_Cholesky *)pc->data; 167 168 PetscFunctionBegin; 169 if (dir->hdr.inplace) { 170 PetscCall(MatSolve(pc->pmat, x, y)); 171 } else { 172 PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y)); 173 } 174 PetscFunctionReturn(0); 175 } 176 177 static PetscErrorCode PCMatApply_Cholesky(PC pc, Mat X, Mat Y) { 178 PC_Cholesky *dir = (PC_Cholesky *)pc->data; 179 180 PetscFunctionBegin; 181 if (dir->hdr.inplace) { 182 PetscCall(MatMatSolve(pc->pmat, X, Y)); 183 } else { 184 PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y)); 185 } 186 PetscFunctionReturn(0); 187 } 188 189 static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc, Vec x, Vec y) { 190 PC_Cholesky *dir = (PC_Cholesky *)pc->data; 191 192 PetscFunctionBegin; 193 if (dir->hdr.inplace) { 194 PetscCall(MatForwardSolve(pc->pmat, x, y)); 195 } else { 196 PetscCall(MatForwardSolve(((PC_Factor *)dir)->fact, x, y)); 197 } 198 PetscFunctionReturn(0); 199 } 200 201 static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc, Vec x, Vec y) { 202 PC_Cholesky *dir = (PC_Cholesky *)pc->data; 203 204 PetscFunctionBegin; 205 if (dir->hdr.inplace) { 206 PetscCall(MatBackwardSolve(pc->pmat, x, y)); 207 } else { 208 PetscCall(MatBackwardSolve(((PC_Factor *)dir)->fact, x, y)); 209 } 210 PetscFunctionReturn(0); 211 } 212 213 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc, Vec x, Vec y) { 214 PC_Cholesky *dir = (PC_Cholesky *)pc->data; 215 216 PetscFunctionBegin; 217 if (dir->hdr.inplace) { 218 PetscCall(MatSolveTranspose(pc->pmat, x, y)); 219 } else { 220 PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y)); 221 } 222 PetscFunctionReturn(0); 223 } 224 225 /*@ 226 PCFactorSetReuseOrdering - When similar matrices are factored, this 227 causes the ordering computed in the first factor to be used for all 228 following factors. 229 230 Logically Collective on pc 231 232 Input Parameters: 233 + pc - the preconditioner context 234 - flag - `PETSC_TRUE` to reuse else `PETSC_FALSE` 235 236 Options Database Key: 237 . -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()` 238 239 Level: intermediate 240 241 .seealso: `PCLU`, `PCCHOLESKY`, `PCFactorSetReuseFill()` 242 @*/ 243 PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag) { 244 PetscFunctionBegin; 245 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 246 PetscValidLogicalCollectiveBool(pc, flag, 2); 247 PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag)); 248 PetscFunctionReturn(0); 249 } 250 251 /*MC 252 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 253 254 Options Database Keys: 255 + -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()` 256 . -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu 257 . -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()` 258 . -pc_factor_fill <fill> - Sets fill amount 259 . -pc_factor_in_place - Activates in-place factorization 260 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 261 262 Level: beginner 263 264 Notes: 265 Not all options work for all matrix formats 266 267 Usually this will compute an "exact" solution in one iteration and does 268 not need a Krylov method (i.e. you can use -ksp_type preonly, or 269 `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method 270 271 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 272 `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`, 273 `PCFactorSetFill()`, `PCFactorSetShiftNonzero()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()` 274 `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()` 275 M*/ 276 277 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc) { 278 PC_Cholesky *dir; 279 280 PetscFunctionBegin; 281 PetscCall(PetscNewLog(pc, &dir)); 282 pc->data = (void *)dir; 283 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY)); 284 285 ((PC_Factor *)dir)->info.fill = 5.0; 286 287 pc->ops->destroy = PCDestroy_Cholesky; 288 pc->ops->reset = PCReset_Cholesky; 289 pc->ops->apply = PCApply_Cholesky; 290 pc->ops->matapply = PCMatApply_Cholesky; 291 pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky; 292 pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky; 293 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 294 pc->ops->setup = PCSetUp_Cholesky; 295 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 296 pc->ops->view = PCView_Factor; 297 pc->ops->applyrichardson = NULL; 298 PetscFunctionReturn(0); 299 } 300