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 227 /* -----------------------------------------------------------------------------------*/ 228 229 /*@ 230 PCFactorSetReuseOrdering - When similar matrices are factored, this 231 causes the ordering computed in the first factor to be used for all 232 following factors. 233 234 Logically Collective on PC 235 236 Input Parameters: 237 + pc - the preconditioner context 238 - flag - PETSC_TRUE to reuse else PETSC_FALSE 239 240 Options Database Key: 241 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 242 243 Level: intermediate 244 245 .seealso: `PCFactorSetReuseFill()` 246 @*/ 247 PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag) { 248 PetscFunctionBegin; 249 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 250 PetscValidLogicalCollectiveBool(pc, flag, 2); 251 PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag)); 252 PetscFunctionReturn(0); 253 } 254 255 /*MC 256 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 257 258 Options Database Keys: 259 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 260 . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu 261 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 262 . -pc_factor_fill <fill> - Sets fill amount 263 . -pc_factor_in_place - Activates in-place factorization 264 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 265 266 Notes: 267 Not all options work for all matrix formats 268 269 Level: beginner 270 271 Notes: 272 Usually this will compute an "exact" solution in one iteration and does 273 not need a Krylov method (i.e. you can use -ksp_type preonly, or 274 KSPSetType(ksp,KSPPREONLY) for the Krylov method 275 276 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 277 `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`, 278 `PCFactorSetFill()`, `PCFactorSetShiftNonzero()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()` 279 `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()` 280 281 M*/ 282 283 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc) { 284 PC_Cholesky *dir; 285 286 PetscFunctionBegin; 287 PetscCall(PetscNewLog(pc, &dir)); 288 pc->data = (void *)dir; 289 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY)); 290 291 ((PC_Factor *)dir)->info.fill = 5.0; 292 293 pc->ops->destroy = PCDestroy_Cholesky; 294 pc->ops->reset = PCReset_Cholesky; 295 pc->ops->apply = PCApply_Cholesky; 296 pc->ops->matapply = PCMatApply_Cholesky; 297 pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky; 298 pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky; 299 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 300 pc->ops->setup = PCSetUp_Cholesky; 301 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 302 pc->ops->view = PCView_Factor; 303 pc->ops->applyrichardson = NULL; 304 PetscFunctionReturn(0); 305 } 306