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