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 8 #include <../src/ksp/pc/impls/factor/lu/lu.h> /*I "petscpc.h" I*/ 9 10 PetscErrorCode PCFactorReorderForNonzeroDiagonal_LU(PC pc, PetscReal z) 11 { 12 PC_LU *lu = (PC_LU *)pc->data; 13 14 PetscFunctionBegin; 15 lu->nonzerosalongdiagonal = PETSC_TRUE; 16 if (z == PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10; 17 else lu->nonzerosalongdiagonaltol = z; 18 PetscFunctionReturn(0); 19 } 20 21 static PetscErrorCode PCSetFromOptions_LU(PC pc, PetscOptionItems *PetscOptionsObject) 22 { 23 PC_LU *lu = (PC_LU *)pc->data; 24 PetscBool flg = PETSC_FALSE; 25 PetscReal tol; 26 27 PetscFunctionBegin; 28 PetscOptionsHeadBegin(PetscOptionsObject, "LU options"); 29 PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject)); 30 31 PetscCall(PetscOptionsName("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", &flg)); 32 if (flg) { 33 tol = PETSC_DECIDE; 34 PetscCall(PetscOptionsReal("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", lu->nonzerosalongdiagonaltol, &tol, NULL)); 35 PetscCall(PCFactorReorderForNonzeroDiagonal(pc, tol)); 36 } 37 PetscOptionsHeadEnd(); 38 PetscFunctionReturn(0); 39 } 40 41 static PetscErrorCode PCSetUp_LU(PC pc) 42 { 43 PC_LU *dir = (PC_LU *)pc->data; 44 MatSolverType stype; 45 MatFactorError err; 46 const char *prefix; 47 48 PetscFunctionBegin; 49 pc->failedreason = PC_NOERROR; 50 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill; 51 52 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 53 PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix)); 54 55 PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure)); 56 if (dir->hdr.inplace) { 57 MatFactorType ftype; 58 59 PetscCall(MatGetFactorType(pc->pmat, &ftype)); 60 if (ftype == MAT_FACTOR_NONE) { 61 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row)); 62 PetscCall(ISDestroy(&dir->col)); 63 /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */ 64 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 65 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col)); 66 if (dir->row) { } 67 PetscCall(MatLUFactor(pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info)); 68 PetscCall(MatFactorGetError(pc->pmat, &err)); 69 if (err) { /* Factor() fails */ 70 pc->failedreason = (PCFailedReason)err; 71 PetscFunctionReturn(0); 72 } 73 } 74 ((PC_Factor *)dir)->fact = pc->pmat; 75 } else { 76 MatInfo info; 77 78 if (!pc->setupcalled) { 79 PetscBool canuseordering; 80 if (!((PC_Factor *)dir)->fact) { PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_LU, &((PC_Factor *)dir)->fact)); } 81 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering)); 82 if (canuseordering) { 83 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 84 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col)); 85 if (dir->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, dir->nonzerosalongdiagonaltol, dir->row, dir->col)); 86 } 87 PetscCall(MatLUFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info)); 88 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info)); 89 dir->hdr.actualfill = info.fill_ratio_needed; 90 } else if (pc->flag != SAME_NONZERO_PATTERN) { 91 PetscBool canuseordering; 92 if (!dir->hdr.reuseordering) { 93 PetscCall(MatDestroy(&((PC_Factor *)dir)->fact)); 94 PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_LU, &((PC_Factor *)dir)->fact)); 95 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering)); 96 if (canuseordering) { 97 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row)); 98 PetscCall(ISDestroy(&dir->col)); 99 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 100 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col)); 101 if (dir->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, dir->nonzerosalongdiagonaltol, dir->row, dir->col)); 102 } 103 } 104 PetscCall(MatLUFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info)); 105 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info)); 106 dir->hdr.actualfill = info.fill_ratio_needed; 107 } else { 108 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 109 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) { 110 PetscCall(MatFactorClearError(((PC_Factor *)dir)->fact)); 111 pc->failedreason = PC_NOERROR; 112 } 113 } 114 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 115 if (err) { /* FactorSymbolic() fails */ 116 pc->failedreason = (PCFailedReason)err; 117 PetscFunctionReturn(0); 118 } 119 120 PetscCall(MatLUFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info)); 121 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 122 if (err) { /* FactorNumeric() fails */ 123 pc->failedreason = (PCFailedReason)err; 124 } 125 } 126 127 PetscCall(PCFactorGetMatSolverType(pc, &stype)); 128 if (!stype) { 129 MatSolverType solverpackage; 130 PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage)); 131 PetscCall(PCFactorSetMatSolverType(pc, solverpackage)); 132 } 133 PetscFunctionReturn(0); 134 } 135 136 static PetscErrorCode PCReset_LU(PC pc) 137 { 138 PC_LU *dir = (PC_LU *)pc->data; 139 140 PetscFunctionBegin; 141 if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact)); 142 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row)); 143 PetscCall(ISDestroy(&dir->col)); 144 PetscFunctionReturn(0); 145 } 146 147 static PetscErrorCode PCDestroy_LU(PC pc) 148 { 149 PC_LU *dir = (PC_LU *)pc->data; 150 151 PetscFunctionBegin; 152 PetscCall(PCReset_LU(pc)); 153 PetscCall(PetscFree(((PC_Factor *)dir)->ordering)); 154 PetscCall(PetscFree(((PC_Factor *)dir)->solvertype)); 155 PetscCall(PCFactorClearComposedFunctions(pc)); 156 PetscCall(PetscFree(pc->data)); 157 PetscFunctionReturn(0); 158 } 159 160 static PetscErrorCode PCApply_LU(PC pc, Vec x, Vec y) 161 { 162 PC_LU *dir = (PC_LU *)pc->data; 163 164 PetscFunctionBegin; 165 if (dir->hdr.inplace) { 166 PetscCall(MatSolve(pc->pmat, x, y)); 167 } else { 168 PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y)); 169 } 170 PetscFunctionReturn(0); 171 } 172 173 static PetscErrorCode PCMatApply_LU(PC pc, Mat X, Mat Y) 174 { 175 PC_LU *dir = (PC_LU *)pc->data; 176 177 PetscFunctionBegin; 178 if (dir->hdr.inplace) { 179 PetscCall(MatMatSolve(pc->pmat, X, Y)); 180 } else { 181 PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y)); 182 } 183 PetscFunctionReturn(0); 184 } 185 186 static PetscErrorCode PCApplyTranspose_LU(PC pc, Vec x, Vec y) 187 { 188 PC_LU *dir = (PC_LU *)pc->data; 189 190 PetscFunctionBegin; 191 if (dir->hdr.inplace) { 192 PetscCall(MatSolveTranspose(pc->pmat, x, y)); 193 } else { 194 PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y)); 195 } 196 PetscFunctionReturn(0); 197 } 198 199 /*MC 200 PCLU - Uses a direct solver, based on LU factorization, as a preconditioner 201 202 Options Database Keys: 203 + -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()` 204 . -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu 205 . -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()` 206 . -pc_factor_fill <fill> - Sets fill amount 207 . -pc_factor_in_place - Activates in-place factorization 208 . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 209 . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase 210 stability of factorization. 211 . -pc_factor_shift_type <shifttype> - Sets shift type or -1 for the default; use '-help' for a list of available types 212 . -pc_factor_shift_amount <shiftamount> - Sets shift amount or -1 for the default 213 . -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the diagonal. 214 . -pc_factor_mat_solver_type <packagename> - use an external package for the solve, see `MatSolverType` for possibilities 215 - -mat_solvertype_optionname - options for a specific solver package, for example -mat_mumps_cntl_1 216 217 Level: beginner 218 219 Notes: 220 Not all options work for all matrix formats 221 222 Run with -help to see additional options for particular matrix formats or factorization algorithms 223 224 Usually this will compute an "exact" solution in one iteration and does 225 not need a Krylov method (i.e. you can use -ksp_type preonly, or 226 `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method 227 228 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MatSolverType`, `MatGetFactor()`, `PCQR`, `PCSVD`, 229 `PCILU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`, 230 `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`, 231 `PCFactorSetPivotInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()` 232 `PCFactorReorderForNonzeroDiagonal()` 233 M*/ 234 235 PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc) 236 { 237 PC_LU *dir; 238 239 PetscFunctionBegin; 240 PetscCall(PetscNew(&dir)); 241 pc->data = (void *)dir; 242 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_LU)); 243 dir->nonzerosalongdiagonal = PETSC_FALSE; 244 245 ((PC_Factor *)dir)->info.fill = 5.0; 246 ((PC_Factor *)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */ 247 ((PC_Factor *)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE; 248 dir->col = NULL; 249 dir->row = NULL; 250 251 pc->ops->reset = PCReset_LU; 252 pc->ops->destroy = PCDestroy_LU; 253 pc->ops->apply = PCApply_LU; 254 pc->ops->matapply = PCMatApply_LU; 255 pc->ops->applytranspose = PCApplyTranspose_LU; 256 pc->ops->setup = PCSetUp_LU; 257 pc->ops->setfromoptions = PCSetFromOptions_LU; 258 pc->ops->view = PCView_Factor; 259 pc->ops->applyrichardson = NULL; 260 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_LU)); 261 PetscFunctionReturn(0); 262 } 263