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,&ltype);
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: }