xref: /petsc/src/ksp/pc/impls/factor/lu/lu.c (revision 5de52c6d2b8d6de382140bd9fae367e1f02f2f13) 
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(PetscOptionItems *PetscOptionsObject,PC pc)
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(PetscOptionsObject,pc));
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   PetscCall(PCGetOptionsPrefix(pc,&prefix));
50   PetscCall(MatSetOptionsPrefix(pc->pmat,prefix));
51   pc->failedreason = PC_NOERROR;
52   if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill;
53 
54   PetscCall(MatSetErrorIfFailure(pc->pmat,pc->erroriffailure));
55   if (dir->hdr.inplace) {
56     MatFactorType ftype;
57 
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       if (dir->row) {
66         PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row));
67         PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col));
68       }
69       PetscCall(MatLUFactor(pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info));
70       PetscCall(MatFactorGetError(pc->pmat,&err));
71       if (err) { /* Factor() fails */
72         pc->failedreason = (PCFailedReason)err;
73         PetscFunctionReturn(0);
74       }
75     }
76     ((PC_Factor*)dir)->fact = pc->pmat;
77   } else {
78     MatInfo info;
79 
80     if (!pc->setupcalled) {
81       PetscBool canuseordering;
82       if (!((PC_Factor*)dir)->fact) {
83         PetscCall(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact));
84         PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact));
85       }
86       PetscCall(MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering));
87       if (canuseordering) {
88         PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
89         PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col));
90         if (dir->nonzerosalongdiagonal) {
91           PetscCall(MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col));
92         }
93         PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row));
94         PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col));
95       }
96       PetscCall(MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info));
97       PetscCall(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info));
98       dir->hdr.actualfill = info.fill_ratio_needed;
99     } else if (pc->flag != SAME_NONZERO_PATTERN) {
100       PetscBool canuseordering;
101       if (!dir->hdr.reuseordering) {
102         PetscCall(MatDestroy(&((PC_Factor*)dir)->fact));
103         PetscCall(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact));
104         PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact));
105         PetscCall(MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering));
106         if (canuseordering) {
107           if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
108           PetscCall(ISDestroy(&dir->col));
109           PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
110           PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col));
111           if (dir->nonzerosalongdiagonal) {
112             PetscCall(MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col));
113           }
114           PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row));
115           PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col));
116         }
117       }
118       PetscCall(MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info));
119       PetscCall(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info));
120       dir->hdr.actualfill = info.fill_ratio_needed;
121     } else {
122       PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err));
123       if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
124         PetscCall(MatFactorClearError(((PC_Factor*)dir)->fact));
125         pc->failedreason = PC_NOERROR;
126       }
127     }
128     PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err));
129     if (err) { /* FactorSymbolic() fails */
130       pc->failedreason = (PCFailedReason)err;
131       PetscFunctionReturn(0);
132     }
133 
134     PetscCall(MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info));
135     PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err));
136     if (err) { /* FactorNumeric() fails */
137       pc->failedreason = (PCFailedReason)err;
138     }
139 
140   }
141 
142   PetscCall(PCFactorGetMatSolverType(pc,&stype));
143   if (!stype) {
144     MatSolverType solverpackage;
145     PetscCall(MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage));
146     PetscCall(PCFactorSetMatSolverType(pc,solverpackage));
147   }
148   PetscFunctionReturn(0);
149 }
150 
151 static PetscErrorCode PCReset_LU(PC pc)
152 {
153   PC_LU          *dir = (PC_LU*)pc->data;
154 
155   PetscFunctionBegin;
156   if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) PetscCall(MatDestroy(&((PC_Factor*)dir)->fact));
157   if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
158   PetscCall(ISDestroy(&dir->col));
159   PetscFunctionReturn(0);
160 }
161 
162 static PetscErrorCode PCDestroy_LU(PC pc)
163 {
164   PC_LU          *dir = (PC_LU*)pc->data;
165 
166   PetscFunctionBegin;
167   PetscCall(PCReset_LU(pc));
168   PetscCall(PetscFree(((PC_Factor*)dir)->ordering));
169   PetscCall(PetscFree(((PC_Factor*)dir)->solvertype));
170   PetscCall(PetscFree(pc->data));
171   PetscFunctionReturn(0);
172 }
173 
174 static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
175 {
176   PC_LU          *dir = (PC_LU*)pc->data;
177 
178   PetscFunctionBegin;
179   if (dir->hdr.inplace) {
180     PetscCall(MatSolve(pc->pmat,x,y));
181   } else {
182     PetscCall(MatSolve(((PC_Factor*)dir)->fact,x,y));
183   }
184   PetscFunctionReturn(0);
185 }
186 
187 static PetscErrorCode PCMatApply_LU(PC pc,Mat X,Mat Y)
188 {
189   PC_LU          *dir = (PC_LU*)pc->data;
190 
191   PetscFunctionBegin;
192   if (dir->hdr.inplace) {
193     PetscCall(MatMatSolve(pc->pmat,X,Y));
194   } else {
195     PetscCall(MatMatSolve(((PC_Factor*)dir)->fact,X,Y));
196   }
197   PetscFunctionReturn(0);
198 }
199 
200 static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
201 {
202   PC_LU          *dir = (PC_LU*)pc->data;
203 
204   PetscFunctionBegin;
205   if (dir->hdr.inplace) {
206     PetscCall(MatSolveTranspose(pc->pmat,x,y));
207   } else {
208     PetscCall(MatSolveTranspose(((PC_Factor*)dir)->fact,x,y));
209   }
210   PetscFunctionReturn(0);
211 }
212 
213 /* -----------------------------------------------------------------------------------*/
214 
215 /*MC
216    PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
217 
218    Options Database Keys:
219 +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
220 .  -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
221 .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
222 .  -pc_factor_fill <fill> - Sets fill amount
223 .  -pc_factor_in_place - Activates in-place factorization
224 .  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
225 .  -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
226                                          stability of factorization.
227 .  -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types
228 .  -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default
229 -   -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the
230         diagonal.
231 
232    Notes:
233     Not all options work for all matrix formats
234           Run with -help to see additional options for particular matrix formats or factorization
235           algorithms
236 
237    Level: beginner
238 
239    Notes:
240     Usually this will compute an "exact" solution in one iteration and does
241           not need a Krylov method (i.e. you can use -ksp_type preonly, or
242           KSPSetType(ksp,KSPPREONLY) for the Krylov method
243 
244 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
245           `PCILU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
246           `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
247           `PCFactorSetPivotInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()`
248           `PCFactorReorderForNonzeroDiagonal()`
249 M*/
250 
251 PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
252 {
253   PC_LU          *dir;
254 
255   PetscFunctionBegin;
256   PetscCall(PetscNewLog(pc,&dir));
257   pc->data = (void*)dir;
258   PetscCall(PCFactorInitialize(pc,MAT_FACTOR_LU));
259   dir->nonzerosalongdiagonal    = PETSC_FALSE;
260 
261   ((PC_Factor*)dir)->info.fill          = 5.0;
262   ((PC_Factor*)dir)->info.dtcol         = 1.e-6;  /* default to pivoting; this is only thing PETSc LU supports */
263   ((PC_Factor*)dir)->info.shifttype     = (PetscReal)MAT_SHIFT_NONE;
264   dir->col                              = NULL;
265   dir->row                              = NULL;
266 
267   pc->ops->reset             = PCReset_LU;
268   pc->ops->destroy           = PCDestroy_LU;
269   pc->ops->apply             = PCApply_LU;
270   pc->ops->matapply          = PCMatApply_LU;
271   pc->ops->applytranspose    = PCApplyTranspose_LU;
272   pc->ops->setup             = PCSetUp_LU;
273   pc->ops->setfromoptions    = PCSetFromOptions_LU;
274   pc->ops->view              = PCView_Factor;
275   pc->ops->applyrichardson   = NULL;
276   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU));
277   PetscFunctionReturn(0);
278 }
279