Actual source code: factimpl.c
petsc-3.7.3 2016-08-01
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->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()");
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: }
24: PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc,PetscReal z)
25: {
26: PC_Factor *ilu = (PC_Factor*)pc->data;
29: ilu->info.zeropivot = z;
30: return(0);
31: }
35: PetscErrorCode PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype)
36: {
37: PC_Factor *dir = (PC_Factor*)pc->data;
40: if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
41: else {
42: dir->info.shifttype = (PetscReal) shifttype;
43: if ((shifttype == MAT_SHIFT_NONZERO || shifttype == MAT_SHIFT_INBLOCKS) && dir->info.shiftamount == 0.0) {
44: dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON; /* set default amount if user has not called PCFactorSetShiftAmount() yet */
45: }
46: }
47: return(0);
48: }
52: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount)
53: {
54: PC_Factor *dir = (PC_Factor*)pc->data;
57: if (shiftamount == (PetscReal) PETSC_DECIDE) dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
58: else dir->info.shiftamount = shiftamount;
59: return(0);
60: }
64: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
65: {
66: PC_Factor *ilu = (PC_Factor*)pc->data;
69: if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
70: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change tolerance after use");
71: }
72: ilu->info.usedt = PETSC_TRUE;
73: ilu->info.dt = dt;
74: ilu->info.dtcol = dtcol;
75: ilu->info.dtcount = dtcount;
76: ilu->info.fill = PETSC_DEFAULT;
77: return(0);
78: }
82: PetscErrorCode PCFactorSetFill_Factor(PC pc,PetscReal fill)
83: {
84: PC_Factor *dir = (PC_Factor*)pc->data;
87: dir->info.fill = fill;
88: return(0);
89: }
93: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc,MatOrderingType ordering)
94: {
95: PC_Factor *dir = (PC_Factor*)pc->data;
97: PetscBool flg;
100: if (!pc->setupcalled) {
101: PetscFree(dir->ordering);
102: PetscStrallocpy(ordering,(char**)&dir->ordering);
103: } else {
104: PetscStrcmp(dir->ordering,ordering,&flg);
105: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change ordering after use");
106: }
107: return(0);
108: }
112: PetscErrorCode PCFactorGetLevels_Factor(PC pc,PetscInt *levels)
113: {
114: PC_Factor *ilu = (PC_Factor*)pc->data;
117: *levels = ilu->info.levels;
118: return(0);
119: }
123: PetscErrorCode PCFactorSetLevels_Factor(PC pc,PetscInt levels)
124: {
125: PC_Factor *ilu = (PC_Factor*)pc->data;
129: if (!pc->setupcalled) ilu->info.levels = levels;
130: else if (ilu->info.levels != levels) {
131: (*pc->ops->reset)(pc); /* remove previous factored matrices */
132: pc->setupcalled = 0; /* force a complete rebuild of preconditioner factored matrices */
133: ilu->info.levels = levels;
134: } else if (ilu->info.usedt) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change levels after use with ILUdt");
135: return(0);
136: }
140: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc,PetscBool flg)
141: {
142: PC_Factor *dir = (PC_Factor*)pc->data;
145: dir->info.diagonal_fill = (PetscReal) flg;
146: return(0);
147: }
151: PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc,PetscBool *flg)
152: {
153: PC_Factor *dir = (PC_Factor*)pc->data;
156: *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
157: return(0);
158: }
160: /* ------------------------------------------------------------------------------------------*/
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: }
175: PetscErrorCode PCFactorGetMatrix_Factor(PC pc,Mat *mat)
176: {
177: PC_Factor *ilu = (PC_Factor*)pc->data;
180: if (!ilu->fact) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
181: *mat = ilu->fact;
182: return(0);
183: }
187: PetscErrorCode PCFactorSetMatSolverPackage_Factor(PC pc,const MatSolverPackage stype)
188: {
190: PC_Factor *lu = (PC_Factor*)pc->data;
193: if (lu->fact) {
194: const MatSolverPackage ltype;
195: PetscBool flg;
196: MatFactorGetSolverPackage(lu->fact,<ype);
197: PetscStrcmp(stype,ltype,&flg);
198: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change solver matrix package after PC has been setup or used");
199: }
200:
201: PetscFree(lu->solvertype);
202: PetscStrallocpy(stype,&lu->solvertype);
203: return(0);
204: }
208: PetscErrorCode PCFactorGetMatSolverPackage_Factor(PC pc,const MatSolverPackage *stype)
209: {
210: PC_Factor *lu = (PC_Factor*)pc->data;
213: *stype = lu->solvertype;
214: return(0);
215: }
219: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol)
220: {
221: PC_Factor *dir = (PC_Factor*)pc->data;
224: 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);
225: dir->info.dtcol = dtcol;
226: return(0);
227: }
231: PetscErrorCode PCSetFromOptions_Factor(PetscOptionItems *PetscOptionsObject,PC pc)
232: {
233: PC_Factor *factor = (PC_Factor*)pc->data;
234: PetscErrorCode ierr;
235: PetscBool flg,set;
236: char tname[256], solvertype[64];
237: PetscFunctionList ordlist;
238: PetscEnum etmp;
239: PetscBool inplace;
242: PCFactorGetUseInPlace(pc,&inplace);
243: PetscOptionsBool("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",inplace,&flg,&set);
244: if (set) {
245: PCFactorSetUseInPlace(pc,flg);
246: }
247: PetscOptionsReal("-pc_factor_fill","Expected non-zeros in factored matrix","PCFactorSetFill",((PC_Factor*)factor)->info.fill,&((PC_Factor*)factor)->info.fill,NULL);
249: PetscOptionsEnum("-pc_factor_shift_type","Type of shift to add to diagonal","PCFactorSetShiftType",MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);
250: if (flg) {
251: PCFactorSetShiftType(pc,(MatFactorShiftType)etmp);
252: }
253: PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,0);
255: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,0);
256: 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);
258: 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);
259: if (set) {
260: PCFactorSetPivotInBlocks(pc,flg);
261: }
263: PetscOptionsBool("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",PETSC_FALSE,&flg,&set);
264: if (set) {
265: PCFactorSetReuseFill(pc,flg);
266: }
267: PetscOptionsBool("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",PETSC_FALSE,&flg,&set);
268: if (set) {
269: PCFactorSetReuseOrdering(pc,flg);
270: }
272: MatGetOrderingList(&ordlist);
273: PetscOptionsFList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in factored matrix","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)factor)->ordering,tname,256,&flg);
274: if (flg) {
275: PCFactorSetMatOrderingType(pc,tname);
276: }
278: /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
279: PetscOptionsString("-pc_factor_mat_solver_package","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,64,&flg);
280: if (flg) {
281: PCFactorSetMatSolverPackage(pc,solvertype);
282: }
283: return(0);
284: }
288: PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer)
289: {
290: PC_Factor *factor = (PC_Factor*)pc->data;
292: PetscBool isstring,iascii;
295: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
296: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
297: if (iascii) {
298: if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
299: if (factor->info.dt > 0) {
300: PetscViewerASCIIPrintf(viewer," drop tolerance %g\n",(double)factor->info.dt);
301: PetscViewerASCIIPrintf(viewer," max nonzeros per row %D\n",factor->info.dtcount);
302: PetscViewerASCIIPrintf(viewer," column permutation tolerance %g\n",(double)factor->info.dtcol);
303: } else if (factor->info.levels == 1) {
304: PetscViewerASCIIPrintf(viewer," %D level of fill\n",(PetscInt)factor->info.levels);
305: } else {
306: PetscViewerASCIIPrintf(viewer," %D levels of fill\n",(PetscInt)factor->info.levels);
307: }
308: }
310: PetscViewerASCIIPrintf(viewer," tolerance for zero pivot %g\n",(double)factor->info.zeropivot);
311: if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */
312: PetscViewerASCIIPrintf(viewer," using %s [%s]\n",MatFactorShiftTypesDetail[(int)factor->info.shifttype],MatFactorShiftTypes[(int)factor->info.shifttype]);
313: }
315: PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",factor->ordering);
317: if (factor->fact) {
318: MatInfo info;
319: MatGetInfo(factor->fact,MAT_LOCAL,&info);
320: PetscViewerASCIIPrintf(viewer," factor fill ratio given %g, needed %g\n",(double)info.fill_ratio_given,(double)info.fill_ratio_needed);
321: PetscViewerASCIIPrintf(viewer," Factored matrix follows:\n");
322: PetscViewerASCIIPushTab(viewer);
323: PetscViewerASCIIPushTab(viewer);
324: PetscViewerASCIIPushTab(viewer);
325: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
326: MatView(factor->fact,viewer);
327: PetscViewerPopFormat(viewer);
328: PetscViewerASCIIPopTab(viewer);
329: PetscViewerASCIIPopTab(viewer);
330: PetscViewerASCIIPopTab(viewer);
331: }
333: } else if (isstring) {
334: MatFactorType t;
335: MatGetFactorType(factor->fact,&t);
336: if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) {
337: PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);
338: }
339: }
340: return(0);
341: }