Actual source code: factimpl.c
2: #include <../src/ksp/pc/impls/factor/factor.h>
4: /* ------------------------------------------------------------------------------------------*/
7: PetscErrorCode PCFactorSetUpMatSolverType_Factor(PC pc)
8: {
9: PC_Factor *icc = (PC_Factor*)pc->data;
13: if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"You can only call this routine after the matrix object has been provided to the solver, for example with KSPSetOperators() or SNESSetJacobian()");
14: if (!pc->setupcalled && !((PC_Factor*)icc)->fact) {
15: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,((PC_Factor*)icc)->factortype,&((PC_Factor*)icc)->fact);
16: }
17: return(0);
18: }
20: PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc,PetscReal z)
21: {
22: PC_Factor *ilu = (PC_Factor*)pc->data;
25: ilu->info.zeropivot = z;
26: return(0);
27: }
29: PetscErrorCode PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype)
30: {
31: PC_Factor *dir = (PC_Factor*)pc->data;
34: if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
35: else {
36: dir->info.shifttype = (PetscReal) shifttype;
37: if ((shifttype == MAT_SHIFT_NONZERO || shifttype == MAT_SHIFT_INBLOCKS) && dir->info.shiftamount == 0.0) {
38: dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON; /* set default amount if user has not called PCFactorSetShiftAmount() yet */
39: }
40: }
41: return(0);
42: }
44: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount)
45: {
46: PC_Factor *dir = (PC_Factor*)pc->data;
49: if (shiftamount == (PetscReal) PETSC_DECIDE) dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
50: else dir->info.shiftamount = shiftamount;
51: return(0);
52: }
54: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
55: {
56: PC_Factor *ilu = (PC_Factor*)pc->data;
59: if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
60: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change tolerance after use");
61: }
62: ilu->info.usedt = PETSC_TRUE;
63: ilu->info.dt = dt;
64: ilu->info.dtcol = dtcol;
65: ilu->info.dtcount = dtcount;
66: ilu->info.fill = PETSC_DEFAULT;
67: return(0);
68: }
70: PetscErrorCode PCFactorSetFill_Factor(PC pc,PetscReal fill)
71: {
72: PC_Factor *dir = (PC_Factor*)pc->data;
75: dir->info.fill = fill;
76: return(0);
77: }
79: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc,MatOrderingType ordering)
80: {
81: PC_Factor *dir = (PC_Factor*)pc->data;
83: PetscBool flg;
86: if (!pc->setupcalled) {
87: PetscFree(dir->ordering);
88: PetscStrallocpy(ordering,(char**)&dir->ordering);
89: } else {
90: PetscStrcmp(dir->ordering,ordering,&flg);
91: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change ordering after use");
92: }
93: return(0);
94: }
96: PetscErrorCode PCFactorGetLevels_Factor(PC pc,PetscInt *levels)
97: {
98: PC_Factor *ilu = (PC_Factor*)pc->data;
101: *levels = ilu->info.levels;
102: return(0);
103: }
105: PetscErrorCode PCFactorGetZeroPivot_Factor(PC pc,PetscReal *pivot)
106: {
107: PC_Factor *ilu = (PC_Factor*)pc->data;
110: *pivot = ilu->info.zeropivot;
111: return(0);
112: }
114: PetscErrorCode PCFactorGetShiftAmount_Factor(PC pc,PetscReal *shift)
115: {
116: PC_Factor *ilu = (PC_Factor*)pc->data;
119: *shift = ilu->info.shiftamount;
120: return(0);
121: }
123: PetscErrorCode PCFactorGetShiftType_Factor(PC pc,MatFactorShiftType *type)
124: {
125: PC_Factor *ilu = (PC_Factor*)pc->data;
128: *type = (MatFactorShiftType) (int) ilu->info.shifttype;
129: return(0);
130: }
132: PetscErrorCode PCFactorSetLevels_Factor(PC pc,PetscInt levels)
133: {
134: PC_Factor *ilu = (PC_Factor*)pc->data;
138: if (!pc->setupcalled) ilu->info.levels = levels;
139: else if (ilu->info.levels != levels) {
140: (*pc->ops->reset)(pc); /* remove previous factored matrices */
141: pc->setupcalled = 0; /* force a complete rebuild of preconditioner factored matrices */
142: ilu->info.levels = levels;
143: } else if (ilu->info.usedt) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change levels after use with ILUdt");
144: return(0);
145: }
147: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc,PetscBool flg)
148: {
149: PC_Factor *dir = (PC_Factor*)pc->data;
152: dir->info.diagonal_fill = (PetscReal) flg;
153: return(0);
154: }
156: PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc,PetscBool *flg)
157: {
158: PC_Factor *dir = (PC_Factor*)pc->data;
161: *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
162: return(0);
163: }
165: /* ------------------------------------------------------------------------------------------*/
167: PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc,PetscBool pivot)
168: {
169: PC_Factor *dir = (PC_Factor*)pc->data;
172: dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
173: return(0);
174: }
176: PetscErrorCode PCFactorGetMatrix_Factor(PC pc,Mat *mat)
177: {
178: PC_Factor *ilu = (PC_Factor*)pc->data;
181: if (!ilu->fact) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
182: *mat = ilu->fact;
183: return(0);
184: }
186: PetscErrorCode PCFactorSetMatSolverType_Factor(PC pc,MatSolverType stype)
187: {
189: PC_Factor *lu = (PC_Factor*)pc->data;
192: if (lu->fact) {
193: MatSolverType ltype;
194: PetscBool flg;
195: MatFactorGetSolverType(lu->fact,<ype);
196: PetscStrcmp(stype,ltype,&flg);
197: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change solver matrix package after PC has been setup or used");
198: }
200: PetscFree(lu->solvertype);
201: PetscStrallocpy(stype,&lu->solvertype);
202: return(0);
203: }
205: PetscErrorCode PCFactorGetMatSolverType_Factor(PC pc,MatSolverType *stype)
206: {
207: PC_Factor *lu = (PC_Factor*)pc->data;
210: *stype = lu->solvertype;
211: return(0);
212: }
214: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol)
215: {
216: PC_Factor *dir = (PC_Factor*)pc->data;
219: if (dtcol < 0.0 || dtcol > 1.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Column pivot tolerance is %g must be between 0 and 1",(double)dtcol);
220: dir->info.dtcol = dtcol;
221: return(0);
222: }
224: PetscErrorCode PCSetFromOptions_Factor(PetscOptionItems *PetscOptionsObject,PC pc)
225: {
226: PC_Factor *factor = (PC_Factor*)pc->data;
227: PetscErrorCode ierr;
228: PetscBool flg,set;
229: char tname[256], solvertype[64];
230: PetscFunctionList ordlist;
231: PetscEnum etmp;
232: PetscBool inplace;
235: PCFactorGetUseInPlace(pc,&inplace);
236: PetscOptionsBool("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",inplace,&flg,&set);
237: if (set) {
238: PCFactorSetUseInPlace(pc,flg);
239: }
240: PetscOptionsReal("-pc_factor_fill","Expected non-zeros in factored matrix","PCFactorSetFill",((PC_Factor*)factor)->info.fill,&((PC_Factor*)factor)->info.fill,NULL);
242: PetscOptionsEnum("-pc_factor_shift_type","Type of shift to add to diagonal","PCFactorSetShiftType",MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);
243: if (flg) {
244: PCFactorSetShiftType(pc,(MatFactorShiftType)etmp);
245: }
246: PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,NULL);
248: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,NULL);
249: 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);
251: PetscOptionsBool("-pc_factor_pivot_in_blocks","Pivot inside matrix dense blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",((PC_Factor*)factor)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
252: if (set) {
253: PCFactorSetPivotInBlocks(pc,flg);
254: }
256: PetscOptionsBool("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",PETSC_FALSE,&flg,&set);
257: if (set) {
258: PCFactorSetReuseFill(pc,flg);
259: }
260: PetscOptionsBool("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",PETSC_FALSE,&flg,&set);
261: if (set) {
262: PCFactorSetReuseOrdering(pc,flg);
263: }
265: MatGetOrderingList(&ordlist);
266: PetscOptionsFList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in factored matrix","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)factor)->ordering,tname,sizeof(tname),&flg);
267: if (flg) {
268: PCFactorSetMatOrderingType(pc,tname);
269: }
271: /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
272: PetscOptionsDeprecated("-pc_factor_mat_solver_package","-pc_factor_mat_solver_type","3.9",NULL);
273: PetscOptionsString("-pc_factor_mat_solver_type","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,sizeof(solvertype),&flg);
274: if (flg) {
275: PCFactorSetMatSolverType(pc,solvertype);
276: }
277: return(0);
278: }
280: PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer)
281: {
282: PC_Factor *factor = (PC_Factor*)pc->data;
283: PetscErrorCode ierr;
284: PetscBool isstring,iascii,flg;
285: MatOrderingType ordering;
288: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
289: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
290: if (iascii) {
291: if (factor->inplace) {
292: PetscViewerASCIIPrintf(viewer," in-place factorization\n");
293: } else {
294: PetscViewerASCIIPrintf(viewer," out-of-place factorization\n");
295: }
297: if (factor->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
298: if (factor->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
299: if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
300: if (factor->info.dt > 0) {
301: PetscViewerASCIIPrintf(viewer," drop tolerance %g\n",(double)factor->info.dt);
302: PetscViewerASCIIPrintf(viewer," max nonzeros per row %D\n",factor->info.dtcount);
303: PetscViewerASCIIPrintf(viewer," column permutation tolerance %g\n",(double)factor->info.dtcol);
304: } else if (factor->info.levels == 1) {
305: PetscViewerASCIIPrintf(viewer," %D level of fill\n",(PetscInt)factor->info.levels);
306: } else {
307: PetscViewerASCIIPrintf(viewer," %D levels of fill\n",(PetscInt)factor->info.levels);
308: }
309: }
311: PetscViewerASCIIPrintf(viewer," tolerance for zero pivot %g\n",(double)factor->info.zeropivot);
312: if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */
313: PetscViewerASCIIPrintf(viewer," using %s [%s]\n",MatFactorShiftTypesDetail[(int)factor->info.shifttype],MatFactorShiftTypes[(int)factor->info.shifttype]);
314: }
316: PetscStrcmp(factor->ordering,MATORDERINGNATURAL_OR_ND,&flg);
317: if (flg) {
318: PetscBool isseqsbaij;
319: PetscObjectTypeCompareAny((PetscObject)pc->pmat,&isseqsbaij,MATSEQSBAIJ,MATSEQBAIJ,NULL);
320: if (isseqsbaij) {
321: ordering = MATORDERINGNATURAL;
322: } else {
323: ordering = MATORDERINGND;
324: }
325: } else {
326: ordering = factor->ordering;
327: }
328: if (!factor->fact) {
329: PetscViewerASCIIPrintf(viewer," matrix ordering: %s (may be overridden during setup)\n",ordering);
330: } else {
331: PetscBool useordering;
332: MatInfo info;
333: MatFactorGetUseOrdering(factor->fact,&useordering);
334: if (!useordering) ordering = "external";
335: PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",ordering);
336: MatGetInfo(factor->fact,MAT_LOCAL,&info);
337: PetscViewerASCIIPrintf(viewer," factor fill ratio given %g, needed %g\n",(double)info.fill_ratio_given,(double)info.fill_ratio_needed);
338: PetscViewerASCIIPrintf(viewer," Factored matrix follows:\n");
339: PetscViewerASCIIPushTab(viewer);
340: PetscViewerASCIIPushTab(viewer);
341: PetscViewerASCIIPushTab(viewer);
342: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
343: MatView(factor->fact,viewer);
344: PetscViewerPopFormat(viewer);
345: PetscViewerASCIIPopTab(viewer);
346: PetscViewerASCIIPopTab(viewer);
347: PetscViewerASCIIPopTab(viewer);
348: }
350: } else if (isstring) {
351: MatFactorType t;
352: MatGetFactorType(factor->fact,&t);
353: if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) {
354: PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);
355: }
356: }
357: return(0);
358: }