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