xref: /petsc/src/ksp/pc/impls/factor/factimpl.c (revision 0700a8246d308f50502909ba325e6169d3ee27eb)
1 #define PETSCKSP_DLL
2 
3 #include "../src/ksp/pc/impls/factor/factor.h"     /*I "petscpc.h"  I*/
4 
5 /* ------------------------------------------------------------------------------------------*/
6 
7 EXTERN_C_BEGIN
8 #undef __FUNCT__
9 #define __FUNCT__ "PCFactorSetZeroPivot_Factor"
10 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetZeroPivot_Factor(PC pc,PetscReal z)
11 {
12   PC_Factor *ilu = (PC_Factor*)pc->data;
13 
14   PetscFunctionBegin;
15   ilu->info.zeropivot = z;
16   PetscFunctionReturn(0);
17 }
18 EXTERN_C_END
19 
20 EXTERN_C_BEGIN
21 #undef __FUNCT__
22 #define __FUNCT__ "PCFactorSetShiftType_Factor"
23 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype)
24 {
25   PC_Factor *dir = (PC_Factor*)pc->data;
26 
27   PetscFunctionBegin;
28   if (shifttype == (MatFactorShiftType)PETSC_DECIDE){
29     dir->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
30   } else {
31     dir->info.shifttype = (PetscReal) shifttype;
32   }
33   PetscFunctionReturn(0);
34 }
35 
36 #undef __FUNCT__
37 #define __FUNCT__ "PCFactorSetShiftAmount_Factor"
38 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount)
39 {
40   PC_Factor *dir = (PC_Factor*)pc->data;
41 
42   PetscFunctionBegin;
43   if (shiftamount == (PetscReal) PETSC_DECIDE){
44     dir->info.shiftamount = 1.e-12;
45   } else {
46     dir->info.shiftamount = shiftamount;
47   }
48   PetscFunctionReturn(0);
49 }
50 EXTERN_C_END
51 
52 EXTERN_C_BEGIN
53 #undef __FUNCT__
54 #define __FUNCT__ "PCFactorSetDropTolerance_Factor"
55 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
56 {
57   PC_Factor         *ilu = (PC_Factor*)pc->data;
58 
59   PetscFunctionBegin;
60   if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
61     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot change tolerance after use");
62   }
63   ilu->info.usedt   = PETSC_TRUE;
64   ilu->info.dt      = dt;
65   ilu->info.dtcol   = dtcol;
66   ilu->info.dtcount = dtcount;
67   ilu->info.fill    = PETSC_DEFAULT;
68   PetscFunctionReturn(0);
69 }
70 EXTERN_C_END
71 
72 EXTERN_C_BEGIN
73 #undef __FUNCT__
74 #define __FUNCT__ "PCFactorSetFill_Factor"
75 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetFill_Factor(PC pc,PetscReal fill)
76 {
77   PC_Factor *dir = (PC_Factor*)pc->data;
78 
79   PetscFunctionBegin;
80   dir->info.fill = fill;
81   PetscFunctionReturn(0);
82 }
83 EXTERN_C_END
84 
85 EXTERN_C_BEGIN
86 #undef __FUNCT__
87 #define __FUNCT__ "PCFactorSetMatOrderingType_Factor"
88 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatOrderingType_Factor(PC pc,const MatOrderingType ordering)
89 {
90   PC_Factor      *dir = (PC_Factor*)pc->data;
91   PetscErrorCode ierr;
92   PetscTruth     flg;
93 
94   PetscFunctionBegin;
95   if (!pc->setupcalled) {
96      ierr = PetscStrfree(dir->ordering);CHKERRQ(ierr);
97      ierr = PetscStrallocpy(ordering,&dir->ordering);CHKERRQ(ierr);
98   } else {
99     ierr = PetscStrcmp(dir->ordering,ordering,&flg);CHKERRQ(ierr);
100     if (!flg) {
101       SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot change ordering after use");
102     }
103   }
104   PetscFunctionReturn(0);
105 }
106 EXTERN_C_END
107 
108 EXTERN_C_BEGIN
109 #undef __FUNCT__
110 #define __FUNCT__ "PCFactorSetLevels_Factor"
111 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetLevels_Factor(PC pc,PetscInt levels)
112 {
113   PC_Factor      *ilu = (PC_Factor*)pc->data;
114 
115   PetscFunctionBegin;
116   if (!pc->setupcalled) {
117     ilu->info.levels = levels;
118   } else if (ilu->info.usedt || ilu->info.levels != levels) {
119     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot change levels after use");
120   }
121   PetscFunctionReturn(0);
122 }
123 EXTERN_C_END
124 
125 EXTERN_C_BEGIN
126 #undef __FUNCT__
127 #define __FUNCT__ "PCFactorSetAllowDiagonalFill_Factor"
128 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetAllowDiagonalFill_Factor(PC pc)
129 {
130   PC_Factor *dir = (PC_Factor*)pc->data;
131 
132   PetscFunctionBegin;
133   dir->info.diagonal_fill = 1;
134   PetscFunctionReturn(0);
135 }
136 EXTERN_C_END
137 
138 
139 /* ------------------------------------------------------------------------------------------*/
140 
141 EXTERN_C_BEGIN
142 #undef __FUNCT__
143 #define __FUNCT__ "PCFactorSetPivotInBlocks_Factor"
144 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetPivotInBlocks_Factor(PC pc,PetscTruth pivot)
145 {
146   PC_Factor *dir = (PC_Factor*)pc->data;
147 
148   PetscFunctionBegin;
149   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
150   PetscFunctionReturn(0);
151 }
152 EXTERN_C_END
153 
154 EXTERN_C_BEGIN
155 #undef __FUNCT__
156 #define __FUNCT__ "PCFactorGetMatrix_Factor"
157 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorGetMatrix_Factor(PC pc,Mat *mat)
158 {
159   PC_Factor *ilu = (PC_Factor*)pc->data;
160 
161   PetscFunctionBegin;
162   if (!ilu->fact) SETERRQ(PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
163   *mat = ilu->fact;
164   PetscFunctionReturn(0);
165 }
166 EXTERN_C_END
167 
168 EXTERN_C_BEGIN
169 #undef __FUNCT__
170 #define __FUNCT__ "PCFactorSetMatSolverPackage_Factor"
171 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatSolverPackage_Factor(PC pc,const MatSolverPackage stype)
172 {
173   PetscErrorCode ierr;
174   PC_Factor      *lu = (PC_Factor*)pc->data;
175 
176   PetscFunctionBegin;
177   if (lu->fact) {
178     const MatSolverPackage ltype;
179     PetscTruth             flg;
180     ierr = MatFactorGetSolverPackage(lu->fact,&ltype);CHKERRQ(ierr);
181     ierr = PetscStrcmp(stype,ltype,&flg);CHKERRQ(ierr);
182     if (!flg) {
183       SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot change solver matrix package after PC has been setup or used");
184     }
185   }
186   ierr = PetscStrfree(lu->solvertype);CHKERRQ(ierr);
187   ierr = PetscStrallocpy(stype,&lu->solvertype);CHKERRQ(ierr);
188   PetscFunctionReturn(0);
189 }
190 EXTERN_C_END
191 
192 EXTERN_C_BEGIN
193 #undef __FUNCT__
194 #define __FUNCT__ "PCFactorGetMatSolverPackage_Factor"
195 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorGetMatSolverPackage_Factor(PC pc,const MatSolverPackage *stype)
196 {
197   PC_Factor      *lu = (PC_Factor*)pc->data;
198 
199   PetscFunctionBegin;
200   *stype = lu->solvertype;
201   PetscFunctionReturn(0);
202 }
203 EXTERN_C_END
204 
205 EXTERN_C_BEGIN
206 #undef __FUNCT__
207 #define __FUNCT__ "PCFactorSetColumnPivot_Factor"
208 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol)
209 {
210   PC_Factor       *dir = (PC_Factor*)pc->data;
211 
212   PetscFunctionBegin;
213   if (dtcol < 0.0 || dtcol > 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Column pivot tolerance is %G must be between 0 and 1",dtcol);
214   dir->info.dtcol = dtcol;
215   PetscFunctionReturn(0);
216 }
217 EXTERN_C_END
218 
219 EXTERN_C_BEGIN
220 #undef __FUNCT__
221 #define __FUNCT__ "PCSetFromOptions_Factor"
222 PetscErrorCode PETSCKSP_DLLEXPORT PCSetFromOptions_Factor(PC pc)
223 {
224   PC_Factor       *factor = (PC_Factor*)pc->data;
225   PetscErrorCode  ierr;
226   PetscTruth      flg = PETSC_FALSE,set;
227   char            tname[256], solvertype[64];
228   PetscFList      ordlist;
229   PetscEnum       etmp;
230 
231   PetscFunctionBegin;
232   if (!MatOrderingRegisterAllCalled) {ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr);}
233   ierr = PetscOptionsTruth("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",flg,&flg,PETSC_NULL);CHKERRQ(ierr);
234   if (flg) {
235     ierr = PCFactorSetUseInPlace(pc);CHKERRQ(ierr);
236   }
237   ierr = PetscOptionsReal("-pc_factor_fill","Expected non-zeros in factored matrix","PCFactorSetFill",((PC_Factor*)factor)->info.fill,&((PC_Factor*)factor)->info.fill,0);CHKERRQ(ierr);
238 
239   ierr = PetscOptionsEnum("-pc_factor_shift_type","Shift added to diagonal","PCFactorSetShiftType",
240                           MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);CHKERRQ(ierr);
241   if (flg) {
242     ((PC_Factor*)factor)->info.shifttype = (PetscReal)(etmp);
243   }
244   ierr = PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,0);CHKERRQ(ierr);
245 
246   ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,0);CHKERRQ(ierr);
247   ierr = PetscOptionsReal("-pc_factor_column_pivot","Column pivot tolerance (used only for some factorization)","PCFactorSetColumnPivot",((PC_Factor*)factor)->info.dtcol,&((PC_Factor*)factor)->info.dtcol,&flg);CHKERRQ(ierr);
248 
249   flg = ((PC_Factor*)factor)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
250   ierr = PetscOptionsTruth("-pc_factor_pivot_in_blocks","Pivot inside matrix dense blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",flg,&flg,&set);CHKERRQ(ierr);
251   if (set) {
252     ierr = PCFactorSetPivotInBlocks(pc,flg);CHKERRQ(ierr);
253   }
254 
255   flg  = PETSC_FALSE;
256   ierr = PetscOptionsTruth("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",flg,&flg,PETSC_NULL);CHKERRQ(ierr);
257   if (flg) {
258     ierr = PCFactorSetReuseFill(pc,PETSC_TRUE);CHKERRQ(ierr);
259   }
260   flg  = PETSC_FALSE;
261   ierr = PetscOptionsTruth("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",flg,&flg,PETSC_NULL);CHKERRQ(ierr);
262   if (flg) {
263     ierr = PCFactorSetReuseOrdering(pc,PETSC_TRUE);CHKERRQ(ierr);
264   }
265 
266   ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr);
267   ierr = PetscOptionsList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in factored matrix","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)factor)->ordering,tname,256,&flg);CHKERRQ(ierr);
268   if (flg) {
269     ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr);
270   }
271 
272   /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
273   ierr = PetscOptionsString("-pc_factor_mat_solver_package","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,64,&flg);CHKERRQ(ierr);
274   if (flg) {
275     ierr = PCFactorSetMatSolverPackage(pc,solvertype);CHKERRQ(ierr);
276   }
277   PetscFunctionReturn(0);
278 }
279 EXTERN_C_END
280 
281 EXTERN_C_BEGIN
282 #include "private/matimpl.h"
283 #undef __FUNCT__
284 #define __FUNCT__ "PCView_Factor"
285 PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer)
286 {
287   PC_Factor       *factor = (PC_Factor*)pc->data;
288   PetscErrorCode  ierr;
289   PetscTruth      isstring,iascii;
290 
291   PetscFunctionBegin;
292   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr);
293   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
294   if (iascii) {
295     if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC){
296       if (factor->info.dt > 0) {
297         ierr = PetscViewerASCIIPrintf(viewer,"  drop tolerance %G\n",factor->info.dt);CHKERRQ(ierr);
298         ierr = PetscViewerASCIIPrintf(viewer,"  max nonzeros per row %D\n",factor->info.dtcount);CHKERRQ(ierr);
299         ierr = PetscViewerASCIIPrintf(viewer,"  column permutation tolerance %G\n",(PetscInt)factor->info.dtcol);CHKERRQ(ierr);
300       } else if (factor->info.levels == 1) {
301         ierr = PetscViewerASCIIPrintf(viewer,"  %D level of fill\n",(PetscInt)factor->info.levels);CHKERRQ(ierr);
302       } else {
303         ierr = PetscViewerASCIIPrintf(viewer,"  %D levels of fill\n",(PetscInt)factor->info.levels);CHKERRQ(ierr);
304       }
305     }
306 
307     ierr = PetscViewerASCIIPrintf(viewer,"  tolerance for zero pivot %G\n",factor->info.zeropivot);CHKERRQ(ierr);
308     if (factor->info.shifttype==(PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
309       ierr = PetscViewerASCIIPrintf(viewer,"  using Manteuffel shift\n");CHKERRQ(ierr);
310     }
311     if (factor->info.shifttype==(PetscReal)MAT_SHIFT_NONZERO) {
312       ierr = PetscViewerASCIIPrintf(viewer,"  using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr);
313     }
314     if (factor->info.shifttype==(PetscReal)MAT_SHIFT_INBLOCKS) {
315       ierr = PetscViewerASCIIPrintf(viewer,"  using diagonal shift on blocks to prevent zero pivot\n");CHKERRQ(ierr);
316     }
317 
318     ierr = PetscViewerASCIIPrintf(viewer,"  matrix ordering: %s\n",factor->ordering);CHKERRQ(ierr);
319 
320     if (factor->fact) {
321       ierr = PetscViewerASCIIPrintf(viewer,"  factor fill ratio given %G, needed %G\n",factor->fact->info.fill_ratio_given,factor->fact->info.fill_ratio_needed);CHKERRQ(ierr);
322       ierr = PetscViewerASCIIPrintf(viewer,"    Factored matrix follows:\n");CHKERRQ(ierr);
323       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
324       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
325       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
326       ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
327       ierr = MatView(factor->fact,viewer);CHKERRQ(ierr);
328       ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
329       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
330       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
331       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
332     }
333 
334   } else if (isstring) {
335     if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC){
336       ierr = PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);CHKERRQ(ierr);CHKERRQ(ierr);
337     }
338   } else {
339     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PC_Factor",((PetscObject)viewer)->type_name);
340   }
341   PetscFunctionReturn(0);
342 }
343 EXTERN_C_END
344