Actual source code: factimpl.c
petsc-3.9.4 2018-09-11
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: }
199:
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,0);
248: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,0);
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,256,&flg);
267: if (flg) {
268: PCFactorSetMatOrderingType(pc,tname);
269: }
271: /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
272: PetscOptionsString("-pc_factor_mat_solver_type","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,64,&flg);
273: if (flg) {
274: PCFactorSetMatSolverType(pc,solvertype);
275: }
276: return(0);
277: }
279: PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer)
280: {
281: PC_Factor *factor = (PC_Factor*)pc->data;
283: PetscBool isstring,iascii;
286: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
287: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
288: if (iascii) {
289: if (factor->inplace) {
290: PetscViewerASCIIPrintf(viewer," in-place factorization\n");
291: } else {
292: PetscViewerASCIIPrintf(viewer," out-of-place factorization\n");
293: }
295: if (factor->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
296: if (factor->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
297: if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
298: if (factor->info.dt > 0) {
299: PetscViewerASCIIPrintf(viewer," drop tolerance %g\n",(double)factor->info.dt);
300: PetscViewerASCIIPrintf(viewer," max nonzeros per row %D\n",factor->info.dtcount);
301: PetscViewerASCIIPrintf(viewer," column permutation tolerance %g\n",(double)factor->info.dtcol);
302: } else if (factor->info.levels == 1) {
303: PetscViewerASCIIPrintf(viewer," %D level of fill\n",(PetscInt)factor->info.levels);
304: } else {
305: PetscViewerASCIIPrintf(viewer," %D levels of fill\n",(PetscInt)factor->info.levels);
306: }
307: }
309: PetscViewerASCIIPrintf(viewer," tolerance for zero pivot %g\n",(double)factor->info.zeropivot);
310: if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */
311: PetscViewerASCIIPrintf(viewer," using %s [%s]\n",MatFactorShiftTypesDetail[(int)factor->info.shifttype],MatFactorShiftTypes[(int)factor->info.shifttype]);
312: }
314: PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",factor->ordering);
316: if (factor->fact) {
317: MatInfo info;
318: MatGetInfo(factor->fact,MAT_LOCAL,&info);
319: PetscViewerASCIIPrintf(viewer," factor fill ratio given %g, needed %g\n",(double)info.fill_ratio_given,(double)info.fill_ratio_needed);
320: PetscViewerASCIIPrintf(viewer," Factored matrix follows:\n");
321: PetscViewerASCIIPushTab(viewer);
322: PetscViewerASCIIPushTab(viewer);
323: PetscViewerASCIIPushTab(viewer);
324: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
325: MatView(factor->fact,viewer);
326: PetscViewerPopFormat(viewer);
327: PetscViewerASCIIPopTab(viewer);
328: PetscViewerASCIIPopTab(viewer);
329: PetscViewerASCIIPopTab(viewer);
330: }
332: } else if (isstring) {
333: MatFactorType t;
334: MatGetFactorType(factor->fact,&t);
335: if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) {
336: PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);
337: }
338: }
339: return(0);
340: }