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