Actual source code: factimpl.c
2: #include <../src/ksp/pc/impls/factor/factor.h>
4: /* ------------------------------------------------------------------------------------------*/
6: PetscErrorCode PCFactorSetUpMatSolverType_Factor(PC pc)
7: {
8: PC_Factor *icc = (PC_Factor*)pc->data;
12: 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()");
13: if (!pc->setupcalled && !((PC_Factor*)icc)->fact) {
14: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,((PC_Factor*)icc)->factortype,&((PC_Factor*)icc)->fact);
15: }
16: return(0);
17: }
19: PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc,PetscReal z)
20: {
21: PC_Factor *ilu = (PC_Factor*)pc->data;
24: ilu->info.zeropivot = z;
25: return(0);
26: }
28: PetscErrorCode PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype)
29: {
30: PC_Factor *dir = (PC_Factor*)pc->data;
33: if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
34: else {
35: dir->info.shifttype = (PetscReal) shifttype;
36: if ((shifttype == MAT_SHIFT_NONZERO || shifttype == MAT_SHIFT_INBLOCKS) && dir->info.shiftamount == 0.0) {
37: dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON; /* set default amount if user has not called PCFactorSetShiftAmount() yet */
38: }
39: }
40: return(0);
41: }
43: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount)
44: {
45: PC_Factor *dir = (PC_Factor*)pc->data;
48: if (shiftamount == (PetscReal) PETSC_DECIDE) dir->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
49: else dir->info.shiftamount = shiftamount;
50: return(0);
51: }
53: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
54: {
55: PC_Factor *ilu = (PC_Factor*)pc->data;
58: if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
59: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change tolerance after use");
60: }
61: ilu->info.usedt = PETSC_TRUE;
62: ilu->info.dt = dt;
63: ilu->info.dtcol = dtcol;
64: ilu->info.dtcount = dtcount;
65: ilu->info.fill = PETSC_DEFAULT;
66: return(0);
67: }
69: PetscErrorCode PCFactorSetFill_Factor(PC pc,PetscReal fill)
70: {
71: PC_Factor *dir = (PC_Factor*)pc->data;
74: dir->info.fill = fill;
75: return(0);
76: }
78: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc,MatOrderingType ordering)
79: {
80: PC_Factor *dir = (PC_Factor*)pc->data;
82: PetscBool flg;
85: if (!pc->setupcalled) {
86: PetscFree(dir->ordering);
87: PetscStrallocpy(ordering,(char**)&dir->ordering);
88: } else {
89: PetscStrcmp(dir->ordering,ordering,&flg);
90: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change ordering after use");
91: }
92: return(0);
93: }
95: PetscErrorCode PCFactorGetLevels_Factor(PC pc,PetscInt *levels)
96: {
97: PC_Factor *ilu = (PC_Factor*)pc->data;
100: *levels = ilu->info.levels;
101: return(0);
102: }
104: PetscErrorCode PCFactorGetZeroPivot_Factor(PC pc,PetscReal *pivot)
105: {
106: PC_Factor *ilu = (PC_Factor*)pc->data;
109: *pivot = ilu->info.zeropivot;
110: return(0);
111: }
113: PetscErrorCode PCFactorGetShiftAmount_Factor(PC pc,PetscReal *shift)
114: {
115: PC_Factor *ilu = (PC_Factor*)pc->data;
118: *shift = ilu->info.shiftamount;
119: return(0);
120: }
122: PetscErrorCode PCFactorGetShiftType_Factor(PC pc,MatFactorShiftType *type)
123: {
124: PC_Factor *ilu = (PC_Factor*)pc->data;
127: *type = (MatFactorShiftType) (int) ilu->info.shifttype;
128: return(0);
129: }
131: PetscErrorCode PCFactorSetLevels_Factor(PC pc,PetscInt levels)
132: {
133: PC_Factor *ilu = (PC_Factor*)pc->data;
137: if (!pc->setupcalled) ilu->info.levels = levels;
138: else if (ilu->info.levels != levels) {
139: (*pc->ops->reset)(pc); /* remove previous factored matrices */
140: pc->setupcalled = 0; /* force a complete rebuild of preconditioner factored matrices */
141: ilu->info.levels = levels;
142: } else if (ilu->info.usedt) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change levels after use with ILUdt");
143: return(0);
144: }
146: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc,PetscBool flg)
147: {
148: PC_Factor *dir = (PC_Factor*)pc->data;
151: dir->info.diagonal_fill = (PetscReal) flg;
152: return(0);
153: }
155: PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc,PetscBool *flg)
156: {
157: PC_Factor *dir = (PC_Factor*)pc->data;
160: *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
161: return(0);
162: }
164: /* ------------------------------------------------------------------------------------------*/
166: PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc,PetscBool pivot)
167: {
168: PC_Factor *dir = (PC_Factor*)pc->data;
171: dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
172: return(0);
173: }
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: }
185: /* allow access to preallocation information */
186: #include <petsc/private/matimpl.h>
188: PetscErrorCode PCFactorSetMatSolverType_Factor(PC pc,MatSolverType stype)
189: {
191: PC_Factor *lu = (PC_Factor*)pc->data;
194: if (lu->fact && lu->fact->assembled) {
195: MatSolverType ltype;
196: PetscBool flg;
197: MatFactorGetSolverType(lu->fact,<ype);
198: PetscStrcmp(stype,ltype,&flg);
199: if (!flg) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change solver matrix package from %s to %s after PC has been setup or used",ltype,stype);
200: }
202: PetscFree(lu->solvertype);
203: PetscStrallocpy(stype,&lu->solvertype);
204: return(0);
205: }
207: PetscErrorCode PCFactorGetMatSolverType_Factor(PC pc,MatSolverType *stype)
208: {
209: PC_Factor *lu = (PC_Factor*)pc->data;
212: *stype = lu->solvertype;
213: return(0);
214: }
216: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol)
217: {
218: PC_Factor *dir = (PC_Factor*)pc->data;
221: 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);
222: dir->info.dtcol = dtcol;
223: return(0);
224: }
226: PetscErrorCode PCSetFromOptions_Factor(PetscOptionItems *PetscOptionsObject,PC pc)
227: {
228: PC_Factor *factor = (PC_Factor*)pc->data;
229: PetscErrorCode ierr;
230: PetscBool flg,set;
231: char tname[256], solvertype[64];
232: PetscFunctionList ordlist;
233: PetscEnum etmp;
234: PetscBool inplace;
237: PCFactorGetUseInPlace(pc,&inplace);
238: PetscOptionsBool("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",inplace,&flg,&set);
239: if (set) {
240: PCFactorSetUseInPlace(pc,flg);
241: }
242: PetscOptionsReal("-pc_factor_fill","Expected non-zeros in factored matrix","PCFactorSetFill",((PC_Factor*)factor)->info.fill,&((PC_Factor*)factor)->info.fill,NULL);
244: PetscOptionsEnum("-pc_factor_shift_type","Type of shift to add to diagonal","PCFactorSetShiftType",MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);
245: if (flg) {
246: PCFactorSetShiftType(pc,(MatFactorShiftType)etmp);
247: }
248: PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,NULL);
250: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,NULL);
251: 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);
253: 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);
254: if (set) {
255: PCFactorSetPivotInBlocks(pc,flg);
256: }
258: PetscOptionsBool("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",PETSC_FALSE,&flg,&set);
259: if (set) {
260: PCFactorSetReuseFill(pc,flg);
261: }
262: PetscOptionsBool("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",PETSC_FALSE,&flg,&set);
263: if (set) {
264: PCFactorSetReuseOrdering(pc,flg);
265: }
267: PetscOptionsDeprecated("-pc_factor_mat_solver_package","-pc_factor_mat_solver_type","3.9",NULL);
268: PetscOptionsString("-pc_factor_mat_solver_type","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,sizeof(solvertype),&flg);
269: if (flg) {
270: PCFactorSetMatSolverType(pc,solvertype);
271: }
272: PCFactorSetDefaultOrdering_Factor(pc);
273: MatGetOrderingList(&ordlist);
274: PetscOptionsFList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in factored matrix","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)factor)->ordering,tname,sizeof(tname),&flg);
275: if (flg) {
276: PCFactorSetMatOrderingType(pc,tname);
277: }
278: return(0);
279: }
281: PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer)
282: {
283: PC_Factor *factor = (PC_Factor*)pc->data;
284: PetscErrorCode ierr;
285: PetscBool isstring,iascii,canuseordering;
286: MatInfo info;
287: MatOrderingType ordering;
290: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
291: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
292: if (iascii) {
293: if (factor->inplace) {
294: PetscViewerASCIIPrintf(viewer," in-place factorization\n");
295: } else {
296: PetscViewerASCIIPrintf(viewer," out-of-place factorization\n");
297: }
299: if (factor->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
300: if (factor->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
301: if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
302: if (factor->info.dt > 0) {
303: PetscViewerASCIIPrintf(viewer," drop tolerance %g\n",(double)factor->info.dt);
304: PetscViewerASCIIPrintf(viewer," max nonzeros per row %D\n",factor->info.dtcount);
305: PetscViewerASCIIPrintf(viewer," column permutation tolerance %g\n",(double)factor->info.dtcol);
306: } else if (factor->info.levels == 1) {
307: PetscViewerASCIIPrintf(viewer," %D level of fill\n",(PetscInt)factor->info.levels);
308: } else {
309: PetscViewerASCIIPrintf(viewer," %D levels of fill\n",(PetscInt)factor->info.levels);
310: }
311: }
313: PetscViewerASCIIPrintf(viewer," tolerance for zero pivot %g\n",(double)factor->info.zeropivot);
314: if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */
315: PetscViewerASCIIPrintf(viewer," using %s [%s]\n",MatFactorShiftTypesDetail[(int)factor->info.shifttype],MatFactorShiftTypes[(int)factor->info.shifttype]);
316: }
318: if (factor->fact) {
319: MatFactorGetCanUseOrdering(factor->fact,&canuseordering);
320: if (!canuseordering) ordering = MATORDERINGEXTERNAL;
321: else ordering = factor->ordering;
322: PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",ordering);
323: if (!factor->fact->assembled) {
324: PetscViewerASCIIPrintf(viewer," matrix solver type: %s\n",factor->fact->solvertype);
325: PetscViewerASCIIPrintf(viewer," matrix not yet factored; no additional information available\n");
326: } else {
327: MatGetInfo(factor->fact,MAT_LOCAL,&info);
328: PetscViewerASCIIPrintf(viewer," factor fill ratio given %g, needed %g\n",(double)info.fill_ratio_given,(double)info.fill_ratio_needed);
329: PetscViewerASCIIPrintf(viewer," Factored matrix follows:\n");
330: PetscViewerASCIIPushTab(viewer);
331: PetscViewerASCIIPushTab(viewer);
332: PetscViewerASCIIPushTab(viewer);
333: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
334: MatView(factor->fact,viewer);
335: PetscViewerPopFormat(viewer);
336: PetscViewerASCIIPopTab(viewer);
337: PetscViewerASCIIPopTab(viewer);
338: PetscViewerASCIIPopTab(viewer);
339: }
340: }
342: } else if (isstring) {
343: MatFactorType t;
344: MatGetFactorType(factor->fact,&t);
345: if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) {
346: PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);
347: }
348: }
349: return(0);
350: }