Actual source code: factimpl.c
petsc-3.6.1 2015-08-06
2: #include <../src/ksp/pc/impls/factor/factor.h> /*I "petscpc.h" I*/
4: /* ------------------------------------------------------------------------------------------*/
9: PetscErrorCode PCFactorSetUpMatSolverPackage_Factor(PC pc)
10: {
11: PC_Factor *icc = (PC_Factor*)pc->data;
15: if (!pc->setupcalled && !((PC_Factor*)icc)->fact) {
16: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,((PC_Factor*)icc)->factortype,&((PC_Factor*)icc)->fact);
17: }
18: return(0);
19: }
23: PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc,PetscReal z)
24: {
25: PC_Factor *ilu = (PC_Factor*)pc->data;
28: ilu->info.zeropivot = z;
29: return(0);
30: }
34: PetscErrorCode PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype)
35: {
36: PC_Factor *dir = (PC_Factor*)pc->data;
39: if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
40: else {
41: dir->info.shifttype = (PetscReal) shifttype;
42: if ((shifttype == MAT_SHIFT_NONZERO || shifttype == MAT_SHIFT_INBLOCKS) && dir->info.shiftamount == 0.0) {
43: dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON; /* set default amount if user has not called PCFactorSetShiftAmount() yet */
44: }
45: }
46: return(0);
47: }
51: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount)
52: {
53: PC_Factor *dir = (PC_Factor*)pc->data;
56: if (shiftamount == (PetscReal) PETSC_DECIDE) dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
57: else dir->info.shiftamount = shiftamount;
58: return(0);
59: }
63: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
64: {
65: PC_Factor *ilu = (PC_Factor*)pc->data;
68: if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
69: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change tolerance after use");
70: }
71: ilu->info.usedt = PETSC_TRUE;
72: ilu->info.dt = dt;
73: ilu->info.dtcol = dtcol;
74: ilu->info.dtcount = dtcount;
75: ilu->info.fill = PETSC_DEFAULT;
76: return(0);
77: }
81: PetscErrorCode PCFactorSetFill_Factor(PC pc,PetscReal fill)
82: {
83: PC_Factor *dir = (PC_Factor*)pc->data;
86: dir->info.fill = fill;
87: return(0);
88: }
92: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc,MatOrderingType ordering)
93: {
94: PC_Factor *dir = (PC_Factor*)pc->data;
96: PetscBool flg;
99: if (!pc->setupcalled) {
100: PetscFree(dir->ordering);
101: PetscStrallocpy(ordering,(char**)&dir->ordering);
102: } else {
103: PetscStrcmp(dir->ordering,ordering,&flg);
104: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change ordering after use");
105: }
106: return(0);
107: }
111: PetscErrorCode PCFactorGetLevels_Factor(PC pc,PetscInt *levels)
112: {
113: PC_Factor *ilu = (PC_Factor*)pc->data;
116: *levels = ilu->info.levels;
117: return(0);
118: }
122: PetscErrorCode PCFactorSetLevels_Factor(PC pc,PetscInt levels)
123: {
124: PC_Factor *ilu = (PC_Factor*)pc->data;
128: if (!pc->setupcalled) ilu->info.levels = levels;
129: else if (ilu->info.levels != levels) {
130: (*pc->ops->reset)(pc); /* remove previous factored matrices */
131: pc->setupcalled = 0; /* force a complete rebuild of preconditioner factored matrices */
132: ilu->info.levels = levels;
133: } else if (ilu->info.usedt) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change levels after use with ILUdt");
134: return(0);
135: }
139: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc,PetscBool flg)
140: {
141: PC_Factor *dir = (PC_Factor*)pc->data;
144: dir->info.diagonal_fill = (PetscReal) flg;
145: return(0);
146: }
150: PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc,PetscBool *flg)
151: {
152: PC_Factor *dir = (PC_Factor*)pc->data;
155: *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
156: return(0);
157: }
159: /* ------------------------------------------------------------------------------------------*/
163: PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc,PetscBool pivot)
164: {
165: PC_Factor *dir = (PC_Factor*)pc->data;
168: dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
169: return(0);
170: }
174: PetscErrorCode PCFactorGetMatrix_Factor(PC pc,Mat *mat)
175: {
176: PC_Factor *ilu = (PC_Factor*)pc->data;
179: if (!ilu->fact) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
180: *mat = ilu->fact;
181: return(0);
182: }
186: PetscErrorCode PCFactorSetMatSolverPackage_Factor(PC pc,const MatSolverPackage stype)
187: {
189: PC_Factor *lu = (PC_Factor*)pc->data;
192: if (lu->fact) {
193: const MatSolverPackage ltype;
194: PetscBool flg;
195: MatFactorGetSolverPackage(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: } else {
199: PetscFree(lu->solvertype);
200: PetscStrallocpy(stype,&lu->solvertype);
201: }
202: return(0);
203: }
207: PetscErrorCode PCFactorGetMatSolverPackage_Factor(PC pc,const MatSolverPackage *stype)
208: {
209: PC_Factor *lu = (PC_Factor*)pc->data;
212: *stype = lu->solvertype;
213: return(0);
214: }
218: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol)
219: {
220: PC_Factor *dir = (PC_Factor*)pc->data;
223: 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);
224: dir->info.dtcol = dtcol;
225: return(0);
226: }
230: PetscErrorCode PCSetFromOptions_Factor(PetscOptions *PetscOptionsObject,PC pc)
231: {
232: PC_Factor *factor = (PC_Factor*)pc->data;
233: PetscErrorCode ierr;
234: PetscBool flg,set;
235: char tname[256], solvertype[64];
236: PetscFunctionList ordlist;
237: PetscEnum etmp;
238: PetscBool inplace;
241: PCFactorGetUseInPlace(pc,&inplace);
242: PetscOptionsBool("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",inplace,&flg,&set);
243: if (set) {
244: PCFactorSetUseInPlace(pc,flg);
245: }
246: PetscOptionsReal("-pc_factor_fill","Expected non-zeros in factored matrix","PCFactorSetFill",((PC_Factor*)factor)->info.fill,&((PC_Factor*)factor)->info.fill,NULL);
248: PetscOptionsEnum("-pc_factor_shift_type","Type of shift to add to diagonal","PCFactorSetShiftType",MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);
249: if (flg) {
250: PCFactorSetShiftType(pc,(MatFactorShiftType)etmp);
251: }
252: PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,0);
254: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,0);
255: 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);
257: 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);
258: if (set) {
259: PCFactorSetPivotInBlocks(pc,flg);
260: }
262: PetscOptionsBool("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",PETSC_FALSE,&flg,&set);
263: if (set) {
264: PCFactorSetReuseFill(pc,flg);
265: }
266: PetscOptionsBool("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",PETSC_FALSE,&flg,&set);
267: if (set) {
268: PCFactorSetReuseOrdering(pc,flg);
269: }
271: MatGetOrderingList(&ordlist);
272: PetscOptionsFList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in factored matrix","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)factor)->ordering,tname,256,&flg);
273: if (flg) {
274: PCFactorSetMatOrderingType(pc,tname);
275: }
277: /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
278: PetscOptionsString("-pc_factor_mat_solver_package","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,64,&flg);
279: if (flg) {
280: PCFactorSetMatSolverPackage(pc,solvertype);
281: }
282: return(0);
283: }
287: PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer)
288: {
289: PC_Factor *factor = (PC_Factor*)pc->data;
291: PetscBool isstring,iascii;
294: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
295: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
296: if (iascii) {
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: }