Actual source code: lu.c
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: */
8: #include <../src/ksp/pc/impls/factor/lu/lu.h>
10: PetscErrorCode PCFactorReorderForNonzeroDiagonal_LU(PC pc, PetscReal z)
11: {
12: PC_LU *lu = (PC_LU *)pc->data;
14: PetscFunctionBegin;
15: lu->nonzerosalongdiagonal = PETSC_TRUE;
16: if (z == (PetscReal)PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10;
17: else lu->nonzerosalongdiagonaltol = z;
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
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;
27: PetscFunctionBegin;
28: PetscOptionsHeadBegin(PetscOptionsObject, "LU options");
29: PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
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(PETSC_SUCCESS);
39: }
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;
48: PetscFunctionBegin;
49: pc->failedreason = PC_NOERROR;
50: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;
52: PetscCall(PCGetOptionsPrefix(pc, &prefix));
53: PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));
55: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
56: if (dir->hdr.inplace) {
57: MatFactorType ftype;
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(PETSC_SUCCESS);
72: }
73: }
74: ((PC_Factor *)dir)->fact = pc->pmat;
75: } else {
76: MatInfo info;
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(PETSC_SUCCESS);
118: }
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: }
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(PETSC_SUCCESS);
134: }
136: static PetscErrorCode PCReset_LU(PC pc)
137: {
138: PC_LU *dir = (PC_LU *)pc->data;
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(PETSC_SUCCESS);
145: }
147: static PetscErrorCode PCDestroy_LU(PC pc)
148: {
149: PC_LU *dir = (PC_LU *)pc->data;
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(PETSC_SUCCESS);
158: }
160: static PetscErrorCode PCApply_LU(PC pc, Vec x, Vec y)
161: {
162: PC_LU *dir = (PC_LU *)pc->data;
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(PETSC_SUCCESS);
171: }
173: static PetscErrorCode PCMatApply_LU(PC pc, Mat X, Mat Y)
174: {
175: PC_LU *dir = (PC_LU *)pc->data;
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(PETSC_SUCCESS);
184: }
186: static PetscErrorCode PCApplyTranspose_LU(PC pc, Vec x, Vec y)
187: {
188: PC_LU *dir = (PC_LU *)pc->data;
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(PETSC_SUCCESS);
197: }
199: /*MC
200: PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
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
217: Level: beginner
219: Notes:
220: Not all options work for all matrix formats
222: Run with -help to see additional options for particular matrix formats or factorization algorithms
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
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*/
235: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
236: {
237: PC_LU *dir;
239: PetscFunctionBegin;
240: PetscCall(PetscNew(&dir));
241: pc->data = (void *)dir;
242: PetscCall(PCFactorInitialize(pc, MAT_FACTOR_LU));
243: dir->nonzerosalongdiagonal = PETSC_FALSE;
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;
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(PETSC_SUCCESS);
262: }