Actual source code: factimpl.c
petsc-3.3-p7 2013-05-11
2: #include <../src/ksp/pc/impls/factor/factor.h> /*I "petscpc.h" I*/
4: /* ------------------------------------------------------------------------------------------*/
7: EXTERN_C_BEGIN
10: PetscErrorCode PCFactorSetUpMatSolverPackage_Factor(PC pc)
11: {
12: PC_Factor *icc = (PC_Factor*)pc->data;
16: if (!pc->setupcalled && !((PC_Factor*)icc)->fact){
17: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,((PC_Factor*)icc)->factortype,& ((PC_Factor*)icc)->fact);
18: }
19: return(0);
20: }
21: EXTERN_C_END
23: EXTERN_C_BEGIN
26: PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc,PetscReal z)
27: {
28: PC_Factor *ilu = (PC_Factor*)pc->data;
31: ilu->info.zeropivot = z;
32: return(0);
33: }
34: EXTERN_C_END
36: EXTERN_C_BEGIN
39: PetscErrorCode PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype)
40: {
41: PC_Factor *dir = (PC_Factor*)pc->data;
44: if (shifttype == (MatFactorShiftType)PETSC_DECIDE){
45: dir->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
46: } else {
47: dir->info.shifttype = (PetscReal) shifttype;
48: if (shifttype == MAT_SHIFT_NONZERO && dir->info.shiftamount == 0.0){
49: dir->info.shiftamount = 1.e-12; /* set default amount if user has not called PCFactorSetShiftAmount() yet */
50: }
51: }
52: return(0);
53: }
57: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount)
58: {
59: PC_Factor *dir = (PC_Factor*)pc->data;
62: if (shiftamount == (PetscReal) PETSC_DECIDE){
63: dir->info.shiftamount = 1.e-12;
64: } else {
65: dir->info.shiftamount = shiftamount;
66: }
67: return(0);
68: }
69: EXTERN_C_END
71: EXTERN_C_BEGIN
74: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
75: {
76: PC_Factor *ilu = (PC_Factor*)pc->data;
79: if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
80: SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change tolerance after use");
81: }
82: ilu->info.usedt = PETSC_TRUE;
83: ilu->info.dt = dt;
84: ilu->info.dtcol = dtcol;
85: ilu->info.dtcount = dtcount;
86: ilu->info.fill = PETSC_DEFAULT;
87: return(0);
88: }
89: EXTERN_C_END
91: EXTERN_C_BEGIN
94: PetscErrorCode PCFactorSetFill_Factor(PC pc,PetscReal fill)
95: {
96: PC_Factor *dir = (PC_Factor*)pc->data;
99: dir->info.fill = fill;
100: return(0);
101: }
102: EXTERN_C_END
104: EXTERN_C_BEGIN
107: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc,const MatOrderingType ordering)
108: {
109: PC_Factor *dir = (PC_Factor*)pc->data;
111: PetscBool flg;
112:
114: if (!pc->setupcalled) {
115: PetscFree(dir->ordering);
116: PetscStrallocpy(ordering,&dir->ordering);
117: } else {
118: PetscStrcmp(dir->ordering,ordering,&flg);
119: if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change ordering after use");
120: }
121: return(0);
122: }
123: EXTERN_C_END
125: EXTERN_C_BEGIN
128: PetscErrorCode PCFactorSetLevels_Factor(PC pc,PetscInt levels)
129: {
130: PC_Factor *ilu = (PC_Factor*)pc->data;
134: if (!pc->setupcalled) {
135: ilu->info.levels = levels;
136: } else if (ilu->info.levels != levels) {
137: (*pc->ops->reset)(pc); /* remove previous factored matrices */
138: pc->setupcalled = 0; /* force a complete rebuild of preconditioner factored matrices */
139: ilu->info.levels = levels;
140: } else if (ilu->info.usedt) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change levels after use with ILUdt");
141: return(0);
142: }
143: EXTERN_C_END
145: EXTERN_C_BEGIN
148: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc)
149: {
150: PC_Factor *dir = (PC_Factor*)pc->data;
153: dir->info.diagonal_fill = 1;
154: return(0);
155: }
156: EXTERN_C_END
159: /* ------------------------------------------------------------------------------------------*/
161: EXTERN_C_BEGIN
164: PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc,PetscBool pivot)
165: {
166: PC_Factor *dir = (PC_Factor*)pc->data;
169: dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
170: return(0);
171: }
172: EXTERN_C_END
176: PetscErrorCode PCFactorGetMatrix_Factor(PC pc,Mat *mat)
177: {
178: PC_Factor *ilu = (PC_Factor*)pc->data;
181: if (!ilu->fact) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
182: *mat = ilu->fact;
183: return(0);
184: }
186: EXTERN_C_BEGIN
189: PetscErrorCode PCFactorSetMatSolverPackage_Factor(PC pc,const MatSolverPackage stype)
190: {
192: PC_Factor *lu = (PC_Factor*)pc->data;
195: if (lu->fact) {
196: const MatSolverPackage ltype;
197: PetscBool flg;
198: MatFactorGetSolverPackage(lu->fact,<ype);
199: PetscStrcmp(stype,ltype,&flg);
200: if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change solver matrix package after PC has been setup or used");
201: } else {
202: PetscFree(lu->solvertype);
203: PetscStrallocpy(stype,&lu->solvertype);
204: }
205: return(0);
206: }
207: EXTERN_C_END
209: EXTERN_C_BEGIN
212: PetscErrorCode PCFactorGetMatSolverPackage_Factor(PC pc,const MatSolverPackage *stype)
213: {
214: PC_Factor *lu = (PC_Factor*)pc->data;
217: *stype = lu->solvertype;
218: return(0);
219: }
220: EXTERN_C_END
222: EXTERN_C_BEGIN
225: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol)
226: {
227: PC_Factor *dir = (PC_Factor*)pc->data;
230: if (dtcol < 0.0 || dtcol > 1.0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column pivot tolerance is %G must be between 0 and 1",dtcol);
231: dir->info.dtcol = dtcol;
232: return(0);
233: }
234: EXTERN_C_END
236: EXTERN_C_BEGIN
239: PetscErrorCode PCSetFromOptions_Factor(PC pc)
240: {
241: PC_Factor *factor = (PC_Factor*)pc->data;
242: PetscErrorCode ierr;
243: PetscBool flg = PETSC_FALSE,set;
244: char tname[256], solvertype[64];
245: PetscFList ordlist;
246: PetscEnum etmp;
249: if (!MatOrderingRegisterAllCalled) {MatOrderingRegisterAll(PETSC_NULL);}
250: PetscOptionsBool("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",flg,&flg,PETSC_NULL);
251: if (flg) {
252: PCFactorSetUseInPlace(pc);
253: }
254: PetscOptionsReal("-pc_factor_fill","Expected non-zeros in factored matrix","PCFactorSetFill",((PC_Factor*)factor)->info.fill,&((PC_Factor*)factor)->info.fill,0);
256: PetscOptionsEnum("-pc_factor_shift_type","Type of shift to add to diagonal","PCFactorSetShiftType",
257: MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);
258: if (flg) {
259: PCFactorSetShiftType(pc,(MatFactorShiftType)etmp);
260: }
261: PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,0);
262:
263: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,0);
264: 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);
266: flg = ((PC_Factor*)factor)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
267: PetscOptionsBool("-pc_factor_pivot_in_blocks","Pivot inside matrix dense blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",flg,&flg,&set);
268: if (set) {
269: PCFactorSetPivotInBlocks(pc,flg);
270: }
271:
272: flg = PETSC_FALSE;
273: PetscOptionsBool("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",flg,&flg,PETSC_NULL);
274: if (flg) {
275: PCFactorSetReuseFill(pc,PETSC_TRUE);
276: }
277: flg = PETSC_FALSE;
278: PetscOptionsBool("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",flg,&flg,PETSC_NULL);
279: if (flg) {
280: PCFactorSetReuseOrdering(pc,PETSC_TRUE);
281: }
283: MatGetOrderingList(&ordlist);
284: PetscOptionsList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in factored matrix","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)factor)->ordering,tname,256,&flg);
285: if (flg) {
286: PCFactorSetMatOrderingType(pc,tname);
287: }
289: /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
290: PetscOptionsString("-pc_factor_mat_solver_package","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,64,&flg);
291: if (flg) {
292: PCFactorSetMatSolverPackage(pc,solvertype);
293: }
294: return(0);
295: }
296: EXTERN_C_END
298: EXTERN_C_BEGIN
301: PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer)
302: {
303: PC_Factor *factor = (PC_Factor*)pc->data;
304: PetscErrorCode ierr;
305: PetscBool isstring,iascii;
308: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
309: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
310: if (iascii) {
311: if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC){
312: if (factor->info.dt > 0) {
313: PetscViewerASCIIPrintf(viewer," drop tolerance %G\n",factor->info.dt);
314: PetscViewerASCIIPrintf(viewer," max nonzeros per row %D\n",factor->info.dtcount);
315: PetscViewerASCIIPrintf(viewer," column permutation tolerance %G\n",(PetscInt)factor->info.dtcol);
316: } else if (factor->info.levels == 1) {
317: PetscViewerASCIIPrintf(viewer," %D level of fill\n",(PetscInt)factor->info.levels);
318: } else {
319: PetscViewerASCIIPrintf(viewer," %D levels of fill\n",(PetscInt)factor->info.levels);
320: }
321: }
322:
323: PetscViewerASCIIPrintf(viewer," tolerance for zero pivot %G\n",factor->info.zeropivot);
324: if (factor->info.shifttype==(PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
325: PetscViewerASCIIPrintf(viewer," using Manteuffel shift\n");
326: }
327: if (factor->info.shifttype==(PetscReal)MAT_SHIFT_NONZERO) {
328: PetscViewerASCIIPrintf(viewer," using diagonal shift to prevent zero pivot\n");
329: }
330: if (factor->info.shifttype==(PetscReal)MAT_SHIFT_INBLOCKS) {
331: PetscViewerASCIIPrintf(viewer," using diagonal shift on blocks to prevent zero pivot\n");
332: }
333:
334: PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",factor->ordering);
335:
336: if (factor->fact) {
337: MatInfo info;
338: MatGetInfo(factor->fact,MAT_LOCAL,&info);
339: PetscViewerASCIIPrintf(viewer," factor fill ratio given %G, needed %G\n",info.fill_ratio_given,info.fill_ratio_needed);
340: PetscViewerASCIIPrintf(viewer," Factored matrix follows:\n");
341: PetscViewerASCIIPushTab(viewer);
342: PetscViewerASCIIPushTab(viewer);
343: PetscViewerASCIIPushTab(viewer);
344: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
345: MatView(factor->fact,viewer);
346: PetscViewerPopFormat(viewer);
347: PetscViewerASCIIPopTab(viewer);
348: PetscViewerASCIIPopTab(viewer);
349: PetscViewerASCIIPopTab(viewer);
350: }
351:
352: } else if (isstring) {
353: MatFactorType t;
354: MatGetFactorType(factor->fact,&t);
355: if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC){
356: PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);
357: }
358: } else {
359: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PC_Factor",((PetscObject)viewer)->type_name);
360: }
361: return(0);
362: }
363: EXTERN_C_END