xref: /petsc/src/ksp/pc/impls/factor/ilu/ilu.c (revision 110e16775a088fcd8eefb98b533461c059f06890) 
1 #define PETSCKSP_DLL
2 
3 /*
4    Defines a ILU factorization preconditioner for any Mat implementation
5 */
6 #include "../src/ksp/pc/impls/factor/ilu/ilu.h"     /*I "petscpc.h"  I*/
7 
8 /* ------------------------------------------------------------------------------------------*/
9 EXTERN_C_BEGIN
10 #undef __FUNCT__
11 #define __FUNCT__ "PCFactorSetReuseFill_ILU"
12 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseFill_ILU(PC pc,PetscTruth flag)
13 {
14   PC_ILU *lu = (PC_ILU*)pc->data;
15 
16   PetscFunctionBegin;
17   lu->reusefill = flag;
18   PetscFunctionReturn(0);
19 }
20 EXTERN_C_END
21 
22 EXTERN_C_BEGIN
23 #undef __FUNCT__
24 #define __FUNCT__ "PCFactorReorderForNonzeroDiagonal_ILU"
25 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)
26 {
27   PC_ILU *ilu = (PC_ILU*)pc->data;
28 
29   PetscFunctionBegin;
30   ilu->nonzerosalongdiagonal = PETSC_TRUE;
31   if (z == PETSC_DECIDE) {
32     ilu->nonzerosalongdiagonaltol = 1.e-10;
33   } else {
34     ilu->nonzerosalongdiagonaltol = z;
35   }
36   PetscFunctionReturn(0);
37 }
38 EXTERN_C_END
39 
40 #undef __FUNCT__
41 #define __FUNCT__ "PCDestroy_ILU_Internal"
42 PetscErrorCode PCDestroy_ILU_Internal(PC pc)
43 {
44   PC_ILU         *ilu = (PC_ILU*)pc->data;
45   PetscErrorCode ierr;
46 
47   PetscFunctionBegin;
48   if (!ilu->inplace && ((PC_Factor*)ilu)->fact) {ierr = MatDestroy(((PC_Factor*)ilu)->fact);CHKERRQ(ierr);}
49   if (ilu->row && ilu->col && ilu->row != ilu->col) {ierr = ISDestroy(ilu->row);CHKERRQ(ierr);}
50   if (ilu->col) {ierr = ISDestroy(ilu->col);CHKERRQ(ierr);}
51   PetscFunctionReturn(0);
52 }
53 
54 EXTERN_C_BEGIN
55 #undef __FUNCT__
56 #define __FUNCT__ "PCFactorSetUseDropTolerance_ILU"
57 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
58 {
59   PC_ILU         *ilu = (PC_ILU*)pc->data;
60   PetscErrorCode ierr;
61 
62   PetscFunctionBegin;
63   if (pc->setupcalled && (!ilu->usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
64     pc->setupcalled   = 0;
65     ierr = PCDestroy_ILU_Internal(pc);CHKERRQ(ierr);
66   }
67   ilu->usedt                      = PETSC_TRUE;
68   ((PC_Factor*)ilu)->info.dt      = dt;
69   ((PC_Factor*)ilu)->info.dtcol   = dtcol;
70   ((PC_Factor*)ilu)->info.dtcount = dtcount;
71   ((PC_Factor*)ilu)->info.fill    = PETSC_DEFAULT;
72   PetscFunctionReturn(0);
73 }
74 EXTERN_C_END
75 
76 EXTERN_C_BEGIN
77 #undef __FUNCT__
78 #define __FUNCT__ "PCFactorSetReuseOrdering_ILU"
79 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering_ILU(PC pc,PetscTruth flag)
80 {
81   PC_ILU *ilu = (PC_ILU*)pc->data;
82 
83   PetscFunctionBegin;
84   ilu->reuseordering = flag;
85   PetscFunctionReturn(0);
86 }
87 EXTERN_C_END
88 
89 EXTERN_C_BEGIN
90 #undef __FUNCT__
91 #define __FUNCT__ "PCFactorSetUseInPlace_ILU"
92 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseInPlace_ILU(PC pc)
93 {
94   PC_ILU *dir = (PC_ILU*)pc->data;
95 
96   PetscFunctionBegin;
97   dir->inplace = PETSC_TRUE;
98   PetscFunctionReturn(0);
99 }
100 EXTERN_C_END
101 
102 #undef __FUNCT__
103 #define __FUNCT__ "PCSetFromOptions_ILU"
104 static PetscErrorCode PCSetFromOptions_ILU(PC pc)
105 {
106   PetscErrorCode ierr;
107   PetscInt       dtmax = 3,itmp;
108   PetscTruth     flg,set;
109   PetscReal      dt[3];
110   char           tname[256], solvertype[64];
111   PC_ILU         *ilu = (PC_ILU*)pc->data;
112   PetscFList     ordlist;
113   PetscReal      tol;
114 
115   PetscFunctionBegin;
116   if (!MatOrderingRegisterAllCalled) {ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr);}
117   ierr = PetscOptionsHead("ILU Options");CHKERRQ(ierr);
118     ierr = PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);CHKERRQ(ierr);
119     if (flg) ((PC_Factor*)ilu)->info.levels = itmp;
120     ierr = PetscOptionsTruth("-pc_factor_in_place","do factorization in place","PCFactorSetUseInPlace",ilu->inplace,&ilu->inplace,PETSC_NULL);CHKERRQ(ierr);
121     flg  = PETSC_FALSE;
122     ierr = PetscOptionsTruth("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",flg,&flg,PETSC_NULL);CHKERRQ(ierr);
123     ((PC_Factor*)ilu)->info.diagonal_fill = (double) flg;
124     ierr = PetscOptionsTruth("-pc_factor_reuse_fill","Reuse fill ratio from previous factorization","PCFactorSetReuseFill",ilu->reusefill,&ilu->reusefill,PETSC_NULL);CHKERRQ(ierr);
125     ierr = PetscOptionsTruth("-pc_factor_reuse_ordering","Reuse previous reordering","PCFactorSetReuseOrdering",ilu->reuseordering,&ilu->reuseordering,PETSC_NULL);CHKERRQ(ierr);
126     flg  = PETSC_FALSE;
127     ierr = PetscOptionsTruth("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",flg,&flg,PETSC_NULL);CHKERRQ(ierr);
128     if (flg) {
129       ierr = PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);CHKERRQ(ierr);
130     }
131     ierr = PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",((PC_Factor*)ilu)->info.shiftnz,&((PC_Factor*)ilu)->info.shiftnz,0);CHKERRQ(ierr);
132     flg = (((PC_Factor*)ilu)->info.shiftpd > 0.0) ? PETSC_TRUE : PETSC_FALSE;
133     ierr = PetscOptionsTruth("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",flg,&flg,PETSC_NULL);CHKERRQ(ierr);
134     ierr = PCFactorSetShiftPd(pc,flg);CHKERRQ(ierr);
135     ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)ilu)->info.zeropivot,&((PC_Factor*)ilu)->info.zeropivot,0);CHKERRQ(ierr);
136 
137     dt[0] = ((PC_Factor*)ilu)->info.dt;
138     dt[1] = ((PC_Factor*)ilu)->info.dtcol;
139     dt[2] = ((PC_Factor*)ilu)->info.dtcount;
140     ierr = PetscOptionsRealArray("-pc_factor_use_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetUseDropTolerance",dt,&dtmax,&flg);CHKERRQ(ierr);
141     if (flg) {
142       ierr = PCFactorSetUseDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);CHKERRQ(ierr);
143     }
144     ierr = PetscOptionsReal("-pc_factor_fill","Expected fill in factorization","PCFactorSetFill",((PC_Factor*)ilu)->info.fill,&((PC_Factor*)ilu)->info.fill,&flg);CHKERRQ(ierr);
145     ierr = PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);CHKERRQ(ierr);
146     if (flg) {
147       tol = PETSC_DECIDE;
148       ierr = PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);CHKERRQ(ierr);
149       ierr = PCFactorReorderForNonzeroDiagonal(pc,tol);CHKERRQ(ierr);
150     }
151 
152     ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr);
153     ierr = PetscOptionsList("-pc_factor_mat_ordering_type","Reorder to reduce nonzeros in ILU","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)ilu)->ordering,tname,256,&flg);CHKERRQ(ierr);
154     if (flg) {
155       ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr);
156     }
157 
158     /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
159     ierr = PetscOptionsString("-pc_factor_mat_solver_package","Specific ILU solver to use","MatGetFactor",((PC_Factor*)ilu)->solvertype,solvertype,64,&flg);CHKERRQ(ierr);
160     if (flg) {
161       ierr = PCFactorSetMatSolverPackage(pc,solvertype);CHKERRQ(ierr);
162     }
163 
164     flg = ((PC_Factor*)ilu)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
165     ierr = PetscOptionsTruth("-pc_factor_pivot_in_blocks","Pivot inside matrix blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",flg,&flg,&set);CHKERRQ(ierr);
166     if (set) {
167       ierr = PCFactorSetPivotInBlocks(pc,flg);CHKERRQ(ierr);
168     }
169     flg  = PETSC_FALSE;
170     ierr = PetscOptionsTruth("-pc_factor_shift_in_blocks","Shift added to diagonal of block","PCFactorSetShiftInBlocks",flg,&flg,PETSC_NULL);CHKERRQ(ierr);
171     if (flg) {
172       ierr = PCFactorSetShiftInBlocks(pc,(PetscReal)PETSC_DECIDE);CHKERRQ(ierr);
173     }
174     ierr = PetscOptionsReal("-pc_factor_shift_in_blocks","Shift added to diagonal of block","PCFactorSetShiftInBlocks",((PC_Factor*)ilu)->info.shiftinblocks,&((PC_Factor*)ilu)->info.shiftinblocks,0);CHKERRQ(ierr);
175 
176   ierr = PetscOptionsTail();CHKERRQ(ierr);
177   PetscFunctionReturn(0);
178 }
179 
180 #undef __FUNCT__
181 #define __FUNCT__ "PCView_ILU"
182 static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
183 {
184   PC_ILU         *ilu = (PC_ILU*)pc->data;
185   PetscErrorCode ierr;
186   PetscTruth     isstring,iascii;
187 
188   PetscFunctionBegin;
189   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr);
190   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
191   if (iascii) {
192     if (ilu->usedt) {
193         ierr = PetscViewerASCIIPrintf(viewer,"  ILU: drop tolerance %G\n",((PC_Factor*)ilu)->info.dt);CHKERRQ(ierr);
194         ierr = PetscViewerASCIIPrintf(viewer,"  ILU: max nonzeros per row %D\n",(PetscInt)((PC_Factor*)ilu)->info.dtcount);CHKERRQ(ierr);
195         ierr = PetscViewerASCIIPrintf(viewer,"  ILU: column permutation tolerance %G\n",((PC_Factor*)ilu)->info.dtcol);CHKERRQ(ierr);
196     } else if (((PC_Factor*)ilu)->info.levels == 1) {
197         ierr = PetscViewerASCIIPrintf(viewer,"  ILU: %D level of fill\n",(PetscInt)((PC_Factor*)ilu)->info.levels);CHKERRQ(ierr);
198     } else {
199         ierr = PetscViewerASCIIPrintf(viewer,"  ILU: %D levels of fill\n",(PetscInt)((PC_Factor*)ilu)->info.levels);CHKERRQ(ierr);
200     }
201     ierr = PetscViewerASCIIPrintf(viewer,"  ILU: factor fill ratio allocated %G\n",((PC_Factor*)ilu)->info.fill);CHKERRQ(ierr);
202     ierr = PetscViewerASCIIPrintf(viewer,"  ILU: tolerance for zero pivot %G\n",((PC_Factor*)ilu)->info.zeropivot);CHKERRQ(ierr);
203     if (((PC_Factor*)ilu)->info.shiftpd) {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: using Manteuffel shift\n");CHKERRQ(ierr);}
204     if (((PC_Factor*)ilu)->info.shiftnz) {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr);}
205     if (((PC_Factor*)ilu)->info.shiftinblocks) {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: using diagonal shift on blocks to prevent zero pivot\n");CHKERRQ(ierr);}
206     if (ilu->inplace) {ierr = PetscViewerASCIIPrintf(viewer,"       in-place factorization\n");CHKERRQ(ierr);}
207     else              {ierr = PetscViewerASCIIPrintf(viewer,"       out-of-place factorization\n");CHKERRQ(ierr);}
208     ierr = PetscViewerASCIIPrintf(viewer,"       matrix ordering: %s\n",((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
209     if (ilu->reusefill)     {ierr = PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorization\n");CHKERRQ(ierr);}
210     if (ilu->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorization\n");CHKERRQ(ierr);}
211     if (((PC_Factor*)ilu)->fact) {
212       ierr = PetscViewerASCIIPrintf(viewer,"  ILU: factor fill ratio needed %G\n",ilu->actualfill);CHKERRQ(ierr);
213       ierr = PetscViewerASCIIPrintf(viewer,"       Factored matrix follows\n");CHKERRQ(ierr);
214       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
215       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
216       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
217       ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
218       ierr = MatView(((PC_Factor*)ilu)->fact,viewer);CHKERRQ(ierr);
219       ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
220       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
221       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
222       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
223     }
224   } else if (isstring) {
225     ierr = PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)((PC_Factor*)ilu)->info.levels,((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);CHKERRQ(ierr);
226   } else {
227     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCILU",((PetscObject)viewer)->type_name);
228   }
229   PetscFunctionReturn(0);
230 }
231 
232 #undef __FUNCT__
233 #define __FUNCT__ "PCSetUp_ILU"
234 static PetscErrorCode PCSetUp_ILU(PC pc)
235 {
236   PetscErrorCode ierr;
237   PC_ILU         *ilu = (PC_ILU*)pc->data;
238   MatInfo        info;
239 
240   PetscFunctionBegin;
241   if (ilu->inplace) {
242     CHKMEMQ;
243     if (!pc->setupcalled) {
244 
245       /* In-place factorization only makes sense with the natural ordering,
246          so we only need to get the ordering once, even if nonzero structure changes */
247       ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
248       if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
249       if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
250     }
251 
252     /* In place ILU only makes sense with fill factor of 1.0 because
253        cannot have levels of fill */
254     ((PC_Factor*)ilu)->info.fill          = 1.0;
255     ((PC_Factor*)ilu)->info.diagonal_fill = 0;
256     ierr = MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
257     CHKMEMQ;
258     ((PC_Factor*)ilu)->fact = pc->pmat;
259   } else if (ilu->usedt) {
260     if (!pc->setupcalled) {
261       ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
262     CHKMEMQ;
263       ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ILUDT,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
264       if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
265       if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
266       ierr = MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
267       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
268     } else if (pc->flag != SAME_NONZERO_PATTERN) {
269     CHKMEMQ;
270       ierr = MatDestroy(((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
271     CHKMEMQ;
272       if (!ilu->reuseordering) {
273         if (ilu->row) {ierr = ISDestroy(ilu->row);CHKERRQ(ierr);}
274         if (ilu->col) {ierr = ISDestroy(ilu->col);CHKERRQ(ierr);}
275         ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
276         if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
277         if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
278       }
279       ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ILUDT,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
280       ierr = MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
281       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
282     } else if (!ilu->reusefill) {
283       ierr = MatDestroy(((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
284       ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ILUDT,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
285       ierr = MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
286       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
287     } else {
288       ierr = MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
289     }
290   } else {
291     if (!pc->setupcalled) {
292       /* first time in so compute reordering and symbolic factorization */
293       ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
294       if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
295       if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
296       /*  Remove zeros along diagonal?     */
297       if (ilu->nonzerosalongdiagonal) {
298         ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr);
299       }
300     CHKMEMQ;
301       ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
302     CHKMEMQ;
303       ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
304       ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
305       ilu->actualfill = info.fill_ratio_needed;
306       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
307     } else if (pc->flag != SAME_NONZERO_PATTERN) {
308       if (!ilu->reuseordering) {
309         /* compute a new ordering for the ILU */
310         ierr = ISDestroy(ilu->row);CHKERRQ(ierr);
311         ierr = ISDestroy(ilu->col);CHKERRQ(ierr);
312         ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
313         if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
314         if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
315         /*  Remove zeros along diagonal?     */
316         if (ilu->nonzerosalongdiagonal) {
317           ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr);
318         }
319       }
320       ierr = MatDestroy(((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
321       ierr = MatGetFactor(pc->pmat,MAT_SOLVER_PETSC,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
322       ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
323       ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
324       ilu->actualfill = info.fill_ratio_needed;
325       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
326     }
327     CHKMEMQ;
328     ierr = MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
329     CHKMEMQ;
330   }
331   PetscFunctionReturn(0);
332 }
333 
334 #undef __FUNCT__
335 #define __FUNCT__ "PCDestroy_ILU"
336 static PetscErrorCode PCDestroy_ILU(PC pc)
337 {
338   PC_ILU         *ilu = (PC_ILU*)pc->data;
339   PetscErrorCode ierr;
340 
341   PetscFunctionBegin;
342   ierr = PCDestroy_ILU_Internal(pc);CHKERRQ(ierr);
343   ierr = PetscStrfree(((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
344   ierr = PetscFree(ilu);CHKERRQ(ierr);
345   PetscFunctionReturn(0);
346 }
347 
348 #undef __FUNCT__
349 #define __FUNCT__ "PCApply_ILU"
350 static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
351 {
352   PC_ILU         *ilu = (PC_ILU*)pc->data;
353   PetscErrorCode ierr;
354 
355   PetscFunctionBegin;
356   ierr = MatSolve(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr);
357   PetscFunctionReturn(0);
358 }
359 
360 #undef __FUNCT__
361 #define __FUNCT__ "PCApplyTranspose_ILU"
362 static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
363 {
364   PC_ILU         *ilu = (PC_ILU*)pc->data;
365   PetscErrorCode ierr;
366 
367   PetscFunctionBegin;
368   ierr = MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr);
369   PetscFunctionReturn(0);
370 }
371 
372 /*MC
373      PCILU - Incomplete factorization preconditioners.
374 
375    Options Database Keys:
376 +  -pc_factor_levels <k> - number of levels of fill for ILU(k)
377 .  -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
378                       its factorization (overwrites original matrix)
379 .  -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
380 .  -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
381 .  -pc_factor_use_drop_tolerance <dt,dtcol,maxrowcount> - use Saad's drop tolerance ILUdt
382 .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
383 .  -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
384                                    this decreases the chance of getting a zero pivot
385 .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
386 .  -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
387                              than 1 the diagonal blocks are factored with partial pivoting (this increases the
388                              stability of the ILU factorization
389 .  -pc_factor_shift_in_blocks - adds a small diagonal to any block if it is singular during ILU factorization
390 .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
391 -  -pc_factor_shift_positive_definite true or false - Activate/Deactivate PCFactorSetShiftPd(); the value
392    is optional with true being the default
393 
394    Level: beginner
395 
396   Concepts: incomplete factorization
397 
398    Notes: Only implemented for some matrix formats. (for parallel use you
399              must use MATMPIROWBS, see MatCreateMPIRowbs(), this supports only ILU(0) and this is not recommended
400              unless you really want a parallel ILU).
401 
402           For BAIJ matrices this implements a point block ILU
403 
404    References:
405    T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
406    self-adjoint elliptic difference equations. SIAM J. Numer. Anal., 5:559--573, 1968.
407 
408    T.A. Oliphant. An implicit numerical method for solving two-dimensional time-dependent dif-
409    fusion problems. Quart. Appl. Math., 19:221--229, 1961.
410 
411    Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST
412       http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical
413       Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
414       Science and Engineering, Kluwer, pp. 167--202.
415 
416 
417 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
418            PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), PCFactorSetUseDropTolerance(),
419            PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
420            PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(),
421            PCFactorSetShiftNonzero(),PCFactorSetShiftPd()
422 
423 M*/
424 
425 EXTERN_C_BEGIN
426 #undef __FUNCT__
427 #define __FUNCT__ "PCCreate_ILU"
428 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ILU(PC pc)
429 {
430   PetscErrorCode ierr;
431   PC_ILU         *ilu;
432 
433   PetscFunctionBegin;
434   ierr = PetscNewLog(pc,PC_ILU,&ilu);CHKERRQ(ierr);
435 
436   ((PC_Factor*)ilu)->fact                    = 0;
437   ierr = MatFactorInfoInitialize(&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
438   ((PC_Factor*)ilu)->info.levels             = 0;
439   ((PC_Factor*)ilu)->info.fill               = 1.0;
440   ilu->col                     = 0;
441   ilu->row                     = 0;
442   ilu->inplace                 = PETSC_FALSE;
443   ierr = PetscStrallocpy(MAT_SOLVER_PETSC,&((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr);
444   ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
445   ilu->reuseordering           = PETSC_FALSE;
446   ilu->usedt                   = PETSC_FALSE;
447   ((PC_Factor*)ilu)->info.dt                 = PETSC_DEFAULT;
448   ((PC_Factor*)ilu)->info.dtcount            = PETSC_DEFAULT;
449   ((PC_Factor*)ilu)->info.dtcol              = PETSC_DEFAULT;
450   ((PC_Factor*)ilu)->info.shiftnz            = 1.e-12;
451   ((PC_Factor*)ilu)->info.shiftpd            = 0.0; /* false */
452   ((PC_Factor*)ilu)->info.zeropivot          = 1.e-12;
453   ((PC_Factor*)ilu)->info.pivotinblocks      = 1.0;
454   ((PC_Factor*)ilu)->info.shiftinblocks      = 1.e-12;
455   ilu->reusefill               = PETSC_FALSE;
456   ((PC_Factor*)ilu)->info.diagonal_fill      = 0;
457   pc->data                     = (void*)ilu;
458 
459   pc->ops->destroy             = PCDestroy_ILU;
460   pc->ops->apply               = PCApply_ILU;
461   pc->ops->applytranspose      = PCApplyTranspose_ILU;
462   pc->ops->setup               = PCSetUp_ILU;
463   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
464   pc->ops->getfactoredmatrix   = PCFactorGetMatrix_Factor;
465   pc->ops->view                = PCView_ILU;
466   pc->ops->applyrichardson     = 0;
467 
468   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",
469                     PCFactorSetZeroPivot_Factor);CHKERRQ(ierr);
470   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Factor",
471                     PCFactorSetShiftNonzero_Factor);CHKERRQ(ierr);
472   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Factor",
473                     PCFactorSetShiftPd_Factor);CHKERRQ(ierr);
474 
475   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",
476                     PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr);
477   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",
478                     PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr);
479   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseDropTolerance_C","PCFactorSetUseDropTolerance_ILU",
480                     PCFactorSetUseDropTolerance_ILU);CHKERRQ(ierr);
481   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",
482                     PCFactorSetFill_Factor);CHKERRQ(ierr);
483   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",
484                     PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr);
485   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_ILU",
486                     PCFactorSetReuseOrdering_ILU);CHKERRQ(ierr);
487   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_ILU",
488                     PCFactorSetReuseFill_ILU);CHKERRQ(ierr);
489   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor",
490                     PCFactorSetLevels_Factor);CHKERRQ(ierr);
491   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_ILU",
492                     PCFactorSetUseInPlace_ILU);CHKERRQ(ierr);
493   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C","PCFactorSetAllowDiagonalFill_Factor",
494                     PCFactorSetAllowDiagonalFill_Factor);CHKERRQ(ierr);
495   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetPivotInBlocks_C","PCFactorSetPivotInBlocks_Factor",
496                     PCFactorSetPivotInBlocks_Factor);CHKERRQ(ierr);
497   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftInBlocks_C","PCFactorSetShiftInBlocks_Factor",
498                     PCFactorSetShiftInBlocks_Factor);CHKERRQ(ierr);
499   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C","PCFactorReorderForNonzeroDiagonal_ILU",
500                     PCFactorReorderForNonzeroDiagonal_ILU);CHKERRQ(ierr);
501   PetscFunctionReturn(0);
502 }
503 EXTERN_C_END
504