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