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