Actual source code: ilu.c

  1: #define PETSCKSP_DLL

  3: /*
  4:    Defines a ILU factorization preconditioner for any Mat implementation
  5: */
 6:  #include private/pcimpl.h
 7:  #include src/ksp/pc/impls/factor/ilu/ilu.h
 8:  #include include/private/matimpl.h

 10: /* ------------------------------------------------------------------------------------------*/
 14: PetscErrorCode  PCFactorSetReuseFill_ILU(PC pc,PetscTruth flag)
 15: {
 16:   PC_ILU *lu;
 17: 
 19:   lu = (PC_ILU*)pc->data;
 20:   lu->reusefill = flag;
 21:   return(0);
 22: }

 28: PetscErrorCode  PCFactorSetZeroPivot_ILU(PC pc,PetscReal z)
 29: {
 30:   PC_ILU *ilu;

 33:   ilu                 = (PC_ILU*)pc->data;
 34:   ilu->info.zeropivot = z;
 35:   return(0);
 36: }

 42: PetscErrorCode  PCFactorSetShiftNonzero_ILU(PC pc,PetscReal shift)
 43: {
 44:   PC_ILU *dir;

 47:   dir = (PC_ILU*)pc->data;
 48:   if (shift == (PetscReal) PETSC_DECIDE) {
 49:     dir->info.shiftnz = 1.e-12;
 50:   } else {
 51:     dir->info.shiftnz = shift;
 52:   }
 53:   return(0);
 54: }

 60: PetscErrorCode  PCFactorSetShiftPd_ILU(PC pc,PetscTruth shift)
 61: {
 62:   PC_ILU *dir;
 63: 
 65:   dir = (PC_ILU*)pc->data;
 66:   if (shift) {
 67:     dir->info.shift_fraction = 0.0;
 68:     dir->info.shiftpd = 1.0;
 69:   } else {
 70:     dir->info.shiftpd = 0.0;
 71:   }
 72:   return(0);
 73: }

 79: PetscErrorCode  PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)
 80: {
 81:   PC_ILU *ilu = (PC_ILU*)pc->data;

 84:   ilu->nonzerosalongdiagonal = PETSC_TRUE;
 85:   if (z == PETSC_DECIDE) {
 86:     ilu->nonzerosalongdiagonaltol = 1.e-10;
 87:   } else {
 88:     ilu->nonzerosalongdiagonaltol = z;
 89:   }
 90:   return(0);
 91: }

 96: PetscErrorCode PCDestroy_ILU_Internal(PC pc)
 97: {
 98:   PC_ILU         *ilu = (PC_ILU*)pc->data;

102:   if (!ilu->inplace && ilu->fact) {MatDestroy(ilu->fact);}
103:   if (ilu->row && ilu->col && ilu->row != ilu->col) {ISDestroy(ilu->row);}
104:   if (ilu->col) {ISDestroy(ilu->col);}
105:   return(0);
106: }

111: PetscErrorCode  PCFactorSetUseDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
112: {
113:   PC_ILU         *ilu;

117:   ilu = (PC_ILU*)pc->data;
118:   if (pc->setupcalled && (!ilu->usedt || ilu->info.dt != dt || ilu->info.dtcol != dtcol || ilu->info.dtcount != dtcount)) {
119:     pc->setupcalled   = 0;
120:     PCDestroy_ILU_Internal(pc);
121:   }
122:   ilu->usedt        = PETSC_TRUE;
123:   ilu->info.dt      = dt;
124:   ilu->info.dtcol   = dtcol;
125:   ilu->info.dtcount = dtcount;
126:   ilu->info.fill    = PETSC_DEFAULT;
127:   return(0);
128: }

134: PetscErrorCode  PCFactorSetFill_ILU(PC pc,PetscReal fill)
135: {
136:   PC_ILU *dir;

139:   dir            = (PC_ILU*)pc->data;
140:   dir->info.fill = fill;
141:   return(0);
142: }

148: PetscErrorCode  PCFactorSetMatOrdering_ILU(PC pc,MatOrderingType ordering)
149: {
150:   PC_ILU         *dir = (PC_ILU*)pc->data;
152:   PetscTruth     flg;
153: 
155:   if (!pc->setupcalled) {
156:      PetscStrfree(dir->ordering);
157:      PetscStrallocpy(ordering,&dir->ordering);
158:   } else {
159:     PetscStrcmp(dir->ordering,ordering,&flg);
160:     if (!flg) {
161:       pc->setupcalled = 0;
162:       PetscStrfree(dir->ordering);
163:       PetscStrallocpy(ordering,&dir->ordering);
164:       /* free the data structures, then create them again */
165:       PCDestroy_ILU_Internal(pc);
166:     }
167:   }
168:   return(0);
169: }

175: PetscErrorCode  PCFactorSetReuseOrdering_ILU(PC pc,PetscTruth flag)
176: {
177:   PC_ILU *ilu;

180:   ilu                = (PC_ILU*)pc->data;
181:   ilu->reuseordering = flag;
182:   return(0);
183: }

189: PetscErrorCode  PCFactorSetLevels_ILU(PC pc,PetscInt levels)
190: {
191:   PC_ILU         *ilu;

195:   ilu = (PC_ILU*)pc->data;

197:   if (!pc->setupcalled) {
198:     ilu->info.levels = levels;
199:   } else if (ilu->usedt || ilu->info.levels != levels) {
200:     ilu->info.levels = levels;
201:     pc->setupcalled  = 0;
202:     ilu->usedt       = PETSC_FALSE;
203:     PCDestroy_ILU_Internal(pc);
204:   }
205:   return(0);
206: }

212: PetscErrorCode  PCFactorSetUseInPlace_ILU(PC pc)
213: {
214:   PC_ILU *dir;

217:   dir          = (PC_ILU*)pc->data;
218:   dir->inplace = PETSC_TRUE;
219:   return(0);
220: }

226: PetscErrorCode  PCFactorSetAllowDiagonalFill_ILU(PC pc)
227: {
228:   PC_ILU *dir;

231:   dir                 = (PC_ILU*)pc->data;
232:   dir->info.diagonal_fill = 1;
233:   return(0);
234: }

239: /*@
240:    PCFactorSetUseDropTolerance - The preconditioner will use an ILU 
241:    based on a drop tolerance.

243:    Collective on PC

245:    Input Parameters:
246: +  pc - the preconditioner context
247: .  dt - the drop tolerance, try from 1.e-10 to .1
248: .  dtcol - tolerance for column pivot, good values [0.1 to 0.01]
249: -  maxrowcount - the max number of nonzeros allowed in a row, best value
250:                  depends on the number of nonzeros in row of original matrix

252:    Options Database Key:
253: .  -pc_factor_use_drop_tolerance <dt,dtcol,maxrowcount> - Sets drop tolerance

255:    Level: intermediate

257:     Notes:
258:       This uses the iludt() code of Saad's SPARSKIT package

260: .keywords: PC, levels, reordering, factorization, incomplete, ILU
261: @*/
262: PetscErrorCode  PCFactorSetUseDropTolerance(PC pc,PetscReal dt,PetscReal dtcol,PetscInt maxrowcount)
263: {
264:   PetscErrorCode ierr,(*f)(PC,PetscReal,PetscReal,PetscInt);

268:   PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetUseDropTolerance_C",(void (**)(void))&f);
269:   if (f) {
270:     (*f)(pc,dt,dtcol,maxrowcount);
271:   }
272:   return(0);
273: }

277: /*@
278:    PCFactorSetLevels - Sets the number of levels of fill to use.

280:    Collective on PC

282:    Input Parameters:
283: +  pc - the preconditioner context
284: -  levels - number of levels of fill

286:    Options Database Key:
287: .  -pc_factor_levels <levels> - Sets fill level

289:    Level: intermediate

291: .keywords: PC, levels, fill, factorization, incomplete, ILU
292: @*/
293: PetscErrorCode  PCFactorSetLevels(PC pc,PetscInt levels)
294: {
295:   PetscErrorCode ierr,(*f)(PC,PetscInt);

299:   if (levels < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"negative levels");
300:   PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetLevels_C",(void (**)(void))&f);
301:   if (f) {
302:     (*f)(pc,levels);
303:   }
304:   return(0);
305: }

309: /*@
310:    PCFactorSetAllowDiagonalFill - Causes all diagonal matrix entries to be 
311:    treated as level 0 fill even if there is no non-zero location.

313:    Collective on PC

315:    Input Parameters:
316: +  pc - the preconditioner context

318:    Options Database Key:
319: .  -pc_factor_diagonal_fill

321:    Notes:
322:    Does not apply with 0 fill.

324:    Level: intermediate

326: .keywords: PC, levels, fill, factorization, incomplete, ILU
327: @*/
328: PetscErrorCode  PCFactorSetAllowDiagonalFill(PC pc)
329: {
330:   PetscErrorCode ierr,(*f)(PC);

334:   PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C",(void (**)(void))&f);
335:   if (f) {
336:     (*f)(pc);
337:   }
338:   return(0);
339: }

341: /* ------------------------------------------------------------------------------------------*/

346: PetscErrorCode  PCFactorSetPivotInBlocks_ILU(PC pc,PetscTruth pivot)
347: {
348:   PC_ILU *dir = (PC_ILU*)pc->data;

351:   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
352:   return(0);
353: }

358: static PetscErrorCode PCSetFromOptions_ILU(PC pc)
359: {
361:   PetscInt       dtmax = 3,itmp;
362:   PetscTruth     flg,set;
363:   PetscReal      dt[3];
364:   char           tname[256];
365:   PC_ILU         *ilu = (PC_ILU*)pc->data;
366:   PetscFList     ordlist;
367:   PetscReal      tol;

370:   MatOrderingRegisterAll(PETSC_NULL);
371:   PetscOptionsHead("ILU Options");
372:     PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)ilu->info.levels,&itmp,&flg);
373:     if (flg) ilu->info.levels = itmp;
374:     PetscOptionsName("-pc_factor_in_place","do factorization in place","PCFactorSetUseInPlace",&flg);
375:     if (flg) ilu->inplace = PETSC_TRUE;
376:     PetscOptionsName("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",&flg);
377:     ilu->info.diagonal_fill = (double) flg;
378:     PetscOptionsName("-pc_factor_reuse_fill","Reuse fill ratio from previous factorization","PCFactorSetReuseFill",&flg);
379:     if (flg) ilu->reusefill = PETSC_TRUE;
380:     PetscOptionsName("-pc_factor_reuse_ordering","Reuse previous reordering","PCFactorSetReuseOrdering",&flg);
381:     if (flg) ilu->reuseordering = PETSC_TRUE;
382:     PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
383:     if (flg) {
384:       PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);
385:     }
386:     PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",ilu->info.shiftnz,&ilu->info.shiftnz,0);
387: 
388:     PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",&flg);
389:     if (flg) {
390:       PetscOptionsInt("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",(PetscInt)ilu->info.shiftpd,&itmp,&flg);
391:       if (flg && !itmp) {
392:         PCFactorSetShiftPd(pc,PETSC_FALSE);
393:       } else {
394:         PCFactorSetShiftPd(pc,PETSC_TRUE);
395:       }
396:     }
397:     PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",ilu->info.zeropivot,&ilu->info.zeropivot,0);

399:     dt[0] = ilu->info.dt;
400:     dt[1] = ilu->info.dtcol;
401:     dt[2] = ilu->info.dtcount;
402:     PetscOptionsRealArray("-pc_factor_use_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetUseDropTolerance",dt,&dtmax,&flg);
403:     if (flg) {
404:       PCFactorSetUseDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
405:     }
406:     PetscOptionsReal("-pc_factor_fill","Expected fill in factorization","PCFactorSetFill",ilu->info.fill,&ilu->info.fill,&flg);
407:     PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
408:     if (flg) {
409:       tol = PETSC_DECIDE;
410:       PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);
411:       PCFactorReorderForNonzeroDiagonal(pc,tol);
412:     }

414:     MatGetOrderingList(&ordlist);
415:     PetscOptionsList("-pc_factor_mat_ordering_type","Reorder to reduce nonzeros in ILU","PCFactorSetMatOrdering",ordlist,ilu->ordering,tname,256,&flg);
416:     if (flg) {
417:       PCFactorSetMatOrdering(pc,tname);
418:     }
419:     flg = ilu->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
420:     PetscOptionsTruth("-pc_factor_pivot_in_blocks","Pivot inside matrix blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",flg,&flg,&set);
421:     if (set) {
422:       PCFactorSetPivotInBlocks(pc,flg);
423:     }
424:   PetscOptionsTail();
425:   return(0);
426: }

430: static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
431: {
432:   PC_ILU         *ilu = (PC_ILU*)pc->data;
434:   PetscTruth     isstring,iascii;

437:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
438:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
439:   if (iascii) {
440:     if (ilu->usedt) {
441:         PetscViewerASCIIPrintf(viewer,"  ILU: drop tolerance %G\n",ilu->info.dt);
442:         PetscViewerASCIIPrintf(viewer,"  ILU: max nonzeros per row %D\n",(PetscInt)ilu->info.dtcount);
443:         PetscViewerASCIIPrintf(viewer,"  ILU: column permutation tolerance %G\n",ilu->info.dtcol);
444:     } else if (ilu->info.levels == 1) {
445:         PetscViewerASCIIPrintf(viewer,"  ILU: %D level of fill\n",(PetscInt)ilu->info.levels);
446:     } else {
447:         PetscViewerASCIIPrintf(viewer,"  ILU: %D levels of fill\n",(PetscInt)ilu->info.levels);
448:     }
449:     PetscViewerASCIIPrintf(viewer,"  ILU: factor fill ratio allocated %G\n",ilu->info.fill);
450:     PetscViewerASCIIPrintf(viewer,"  ILU: tolerance for zero pivot %G\n",ilu->info.zeropivot);
451:     if (ilu->info.shiftpd) {PetscViewerASCIIPrintf(viewer,"  ILU: using Manteuffel shift\n");}
452:     if (ilu->inplace) {PetscViewerASCIIPrintf(viewer,"       in-place factorization\n");}
453:     else              {PetscViewerASCIIPrintf(viewer,"       out-of-place factorization\n");}
454:     PetscViewerASCIIPrintf(viewer,"       matrix ordering: %s\n",ilu->ordering);
455:     if (ilu->reusefill)     {PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorization\n");}
456:     if (ilu->reuseordering) {PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorization\n");}
457:     if (ilu->fact) {
458:       PetscViewerASCIIPrintf(viewer,"  ILU: factor fill ratio needed %G\n",ilu->actualfill);
459:       PetscViewerASCIIPrintf(viewer,"       Factored matrix follows\n");
460:       PetscViewerASCIIPushTab(viewer);
461:       PetscViewerASCIIPushTab(viewer);
462:       PetscViewerASCIIPushTab(viewer);
463:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
464:       MatView(ilu->fact,viewer);
465:       PetscViewerPopFormat(viewer);
466:       PetscViewerASCIIPopTab(viewer);
467:       PetscViewerASCIIPopTab(viewer);
468:       PetscViewerASCIIPopTab(viewer);
469:     }
470:   } else if (isstring) {
471:     PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)ilu->info.levels,ilu->ordering);
472:   } else {
473:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCILU",((PetscObject)viewer)->type_name);
474:   }
475:   return(0);
476: }

480: static PetscErrorCode PCSetUp_ILU(PC pc)
481: {
483:   PC_ILU         *ilu = (PC_ILU*)pc->data;
484:   MatInfo        info;

487:   if (ilu->inplace) {
488:     if (!pc->setupcalled) {

490:       /* In-place factorization only makes sense with the natural ordering,
491:          so we only need to get the ordering once, even if nonzero structure changes */
492:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
493:       if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
494:       if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
495:     }

497:     /* In place ILU only makes sense with fill factor of 1.0 because 
498:        cannot have levels of fill */
499:     ilu->info.fill          = 1.0;
500:     ilu->info.diagonal_fill = 0;
501:     MatILUFactor(pc->pmat,ilu->row,ilu->col,&ilu->info);
502:     ilu->fact = pc->pmat;
503:   } else if (ilu->usedt) {
504:     if (!pc->setupcalled) {
505:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
506:       if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
507:       if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
508:       MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
509:       PetscLogObjectParent(pc,ilu->fact);
510:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
511:       MatDestroy(ilu->fact);
512:       if (!ilu->reuseordering) {
513:         if (ilu->row) {ISDestroy(ilu->row);}
514:         if (ilu->col) {ISDestroy(ilu->col);}
515:         MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
516:         if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
517:         if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
518:       }
519:       MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
520:       PetscLogObjectParent(pc,ilu->fact);
521:     } else if (!ilu->reusefill) {
522:       MatDestroy(ilu->fact);
523:       MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
524:       PetscLogObjectParent(pc,ilu->fact);
525:     } else {
526:       MatLUFactorNumeric(pc->pmat,&ilu->info,&ilu->fact);
527:     }
528:   } else {
529:     if (!pc->setupcalled) {
530:       /* first time in so compute reordering and symbolic factorization */
531:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
532:       if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
533:       if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
534:       /*  Remove zeros along diagonal?     */
535:       if (ilu->nonzerosalongdiagonal) {
536:         MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
537:       }
538:       MatILUFactorSymbolic(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
539:       MatGetInfo(ilu->fact,MAT_LOCAL,&info);
540:       ilu->actualfill = info.fill_ratio_needed;
541:       PetscLogObjectParent(pc,ilu->fact);
542:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
543:       if (!ilu->reuseordering) {
544:         /* compute a new ordering for the ILU */
545:         ISDestroy(ilu->row);
546:         ISDestroy(ilu->col);
547:         MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
548:         if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
549:         if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
550:         /*  Remove zeros along diagonal?     */
551:         if (ilu->nonzerosalongdiagonal) {
552:           MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
553:         }
554:       }
555:       MatDestroy(ilu->fact);
556:       MatILUFactorSymbolic(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
557:       MatGetInfo(ilu->fact,MAT_LOCAL,&info);
558:       ilu->actualfill = info.fill_ratio_needed;
559:       PetscLogObjectParent(pc,ilu->fact);
560:     }
561:     MatLUFactorNumeric(pc->pmat,&ilu->info,&ilu->fact);
562:   }
563:   return(0);
564: }

568: static PetscErrorCode PCDestroy_ILU(PC pc)
569: {
570:   PC_ILU         *ilu = (PC_ILU*)pc->data;

574:   PCDestroy_ILU_Internal(pc);
575:   PetscStrfree(ilu->ordering);
576:   PetscFree(ilu);
577:   return(0);
578: }

582: static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
583: {
584:   PC_ILU         *ilu = (PC_ILU*)pc->data;

588:   MatSolve(ilu->fact,x,y);
589:   return(0);
590: }

594: static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
595: {
596:   PC_ILU         *ilu = (PC_ILU*)pc->data;

600:   MatSolveTranspose(ilu->fact,x,y);
601:   return(0);
602: }

606: static PetscErrorCode PCGetFactoredMatrix_ILU(PC pc,Mat *mat)
607: {
608:   PC_ILU *ilu = (PC_ILU*)pc->data;

611:   if (!ilu->fact) SETERRQ(PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
612:   *mat = ilu->fact;
613:   return(0);
614: }

616: /*MC
617:      PCILU - Incomplete factorization preconditioners.

619:    Options Database Keys:
620: +  -pc_factor_levels <k> - number of levels of fill for ILU(k)
621: .  -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
622:                       its factorization (overwrites original matrix)
623: .  -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
624: .  -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
625: .  -pc_factor_use_drop_tolerance <dt,dtcol,maxrowcount> - use Saad's drop tolerance ILUdt
626: .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
627: .  -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
628:                                    this decreases the chance of getting a zero pivot
629: .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
630: .  -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
631:                              than 1 the diagonal blocks are factored with partial pivoting (this increases the 
632:                              stability of the ILU factorization
633: .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
634: -  -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
635:    is optional with PETSC_TRUE being the default

637:    Level: beginner

639:   Concepts: incomplete factorization

641:    Notes: Only implemented for some matrix formats. Not implemented in parallel

643:           For BAIJ matrices this implements a point block ILU

645: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
646:            PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), PCFactorSetUseDropTolerance(),
647:            PCFactorSetFill(), PCFactorSetMatOrdering(), PCFactorSetReuseOrdering(),
648:            PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(),
649:            PCFactorSetShiftNonzero(),PCFactorSetShiftPd()

651: M*/

656: PetscErrorCode  PCCreate_ILU(PC pc)
657: {
659:   PC_ILU         *ilu;

662:   PetscNew(PC_ILU,&ilu);
663:   PetscLogObjectMemory(pc,sizeof(PC_ILU));

665:   ilu->fact                    = 0;
666:   MatFactorInfoInitialize(&ilu->info);
667:   ilu->info.levels             = 0;
668:   ilu->info.fill               = 1.0;
669:   ilu->col                     = 0;
670:   ilu->row                     = 0;
671:   ilu->inplace                 = PETSC_FALSE;
672:   PetscStrallocpy(MATORDERING_NATURAL,&ilu->ordering);
673:   ilu->reuseordering           = PETSC_FALSE;
674:   ilu->usedt                   = PETSC_FALSE;
675:   ilu->info.dt                 = PETSC_DEFAULT;
676:   ilu->info.dtcount            = PETSC_DEFAULT;
677:   ilu->info.dtcol              = PETSC_DEFAULT;
678:   ilu->info.shiftnz            = 0.0;
679:   ilu->info.shiftpd            = 0.0; /* false */
680:   ilu->info.shift_fraction     = 0.0;
681:   ilu->info.zeropivot          = 1.e-12;
682:   ilu->info.pivotinblocks      = 1.0;
683:   ilu->reusefill               = PETSC_FALSE;
684:   ilu->info.diagonal_fill      = 0;
685:   pc->data                     = (void*)ilu;

687:   pc->ops->destroy             = PCDestroy_ILU;
688:   pc->ops->apply               = PCApply_ILU;
689:   pc->ops->applytranspose      = PCApplyTranspose_ILU;
690:   pc->ops->setup               = PCSetUp_ILU;
691:   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
692:   pc->ops->getfactoredmatrix   = PCGetFactoredMatrix_ILU;
693:   pc->ops->view                = PCView_ILU;
694:   pc->ops->applyrichardson     = 0;

696:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_ILU",
697:                     PCFactorSetZeroPivot_ILU);
698:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_ILU",
699:                     PCFactorSetShiftNonzero_ILU);
700:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_ILU",
701:                     PCFactorSetShiftPd_ILU);

703:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseDropTolerance_C","PCFactorSetUseDropTolerance_ILU",
704:                     PCFactorSetUseDropTolerance_ILU);
705:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_ILU",
706:                     PCFactorSetFill_ILU);
707:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrdering_C","PCFactorSetMatOrdering_ILU",
708:                     PCFactorSetMatOrdering_ILU);
709:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_ILU",
710:                     PCFactorSetReuseOrdering_ILU);
711:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_ILU",
712:                     PCFactorSetReuseFill_ILU);
713:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_ILU",
714:                     PCFactorSetLevels_ILU);
715:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_ILU",
716:                     PCFactorSetUseInPlace_ILU);
717:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C","PCFactorSetAllowDiagonalFill_ILU",
718:                     PCFactorSetAllowDiagonalFill_ILU);
719:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetPivotInBlocks_C","PCFactorSetPivotInBlocks_ILU",
720:                     PCFactorSetPivotInBlocks_ILU);
721:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C","PCFactorReorderForNonzeroDiagonal_ILU",
722:                     PCFactorReorderForNonzeroDiagonal_ILU);
723:   return(0);
724: }