Actual source code: ilu.c

  1: /*$Id: ilu.c,v 1.172 2001/08/07 21:30:27 bsmith Exp $*/
  2: /*
  3:    Defines a ILU factorization preconditioner for any Mat implementation
  4: */
 5:  #include src/sles/pc/pcimpl.h
 6:  #include src/sles/pc/impls/ilu/ilu.h
 7:  #include src/mat/matimpl.h

  9: /* ------------------------------------------------------------------------------------------*/

 11: EXTERN_C_BEGIN
 12: #undef __FUNCT__  
 14: int PCILUSetZeroPivot_ILU(PC pc,PetscReal z)
 15: {
 16:   PC_ILU *lu;

 19:   lu                 = (PC_ILU*)pc->data;
 20:   lu->info.zeropivot = z;
 21:   return(0);
 22: }
 23: EXTERN_C_END

 25: #undef __FUNCT__  
 27: int PCDestroy_ILU_Internal(PC pc)
 28: {
 29:   PC_ILU *ilu = (PC_ILU*)pc->data;
 30:   int    ierr;

 33:   if (!ilu->inplace && ilu->fact) {MatDestroy(ilu->fact);}
 34:   if (ilu->row && ilu->col && ilu->row != ilu->col) {ISDestroy(ilu->row);}
 35:   if (ilu->col) {ISDestroy(ilu->col);}
 36:   return(0);
 37: }

 39: EXTERN_C_BEGIN
 40: #undef __FUNCT__  
 42: int PCILUSetDamping_ILU(PC pc,PetscReal damping)
 43: {
 44:   PC_ILU *dir;

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

 57: EXTERN_C_BEGIN
 58: #undef __FUNCT__  
 60: int PCILUSetUseDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,int dtcount)
 61: {
 62:   PC_ILU *ilu;
 63:   int    ierr;

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

 68:   if (!pc->setupcalled) {
 69:     ilu->usedt        = PETSC_TRUE;
 70:     ilu->info.dt      = dt;
 71:     ilu->info.dtcol   = dtcol;
 72:     ilu->info.dtcount = dtcount;
 73:     ilu->info.fill    = PETSC_DEFAULT;
 74:   } else if ((ilu->usedt == PETSC_FALSE)
 75:     || ilu->info.dt != dt || ilu->info.dtcol != dtcol || ilu->info.dtcount != dtcount) {
 76:     ilu->usedt        = PETSC_TRUE;
 77:     ilu->info.dt      = dt;
 78:     ilu->info.dtcol   = dtcol;
 79:     ilu->info.dtcount = dtcount;
 80:     ilu->info.fill    = PETSC_DEFAULT;
 81:     pc->setupcalled   = 0;
 82:     PCDestroy_ILU_Internal(pc);
 83:   }
 84:   return(0);
 85: }
 86: EXTERN_C_END

 88: EXTERN_C_BEGIN
 89: #undef __FUNCT__  
 91: int PCILUSetFill_ILU(PC pc,PetscReal fill)
 92: {
 93:   PC_ILU *dir;

 96:   dir            = (PC_ILU*)pc->data;
 97:   dir->info.fill = fill;
 98:   return(0);
 99: }
100: EXTERN_C_END

102: EXTERN_C_BEGIN
103: #undef __FUNCT__  
105: int PCILUSetMatOrdering_ILU(PC pc,MatOrderingType ordering)
106: {
107:   PC_ILU *dir = (PC_ILU*)pc->data;
108:   int    ierr;
109: 
111:   if (!pc->setupcalled) {
112:      PetscStrfree(dir->ordering);
113:      PetscStrallocpy(ordering,&dir->ordering);
114:   } else if (dir->ordering != ordering) {
115:      pc->setupcalled = 0;
116:      PetscStrfree(dir->ordering);
117:      PetscStrallocpy(ordering,&dir->ordering);
118:      /* free the data structures, then create them again */
119:      PCDestroy_ILU_Internal(pc);
120:   }
121:   return(0);
122: }
123: EXTERN_C_END

125: EXTERN_C_BEGIN
126: #undef __FUNCT__  
128: int PCILUSetReuseOrdering_ILU(PC pc,PetscTruth flag)
129: {
130:   PC_ILU *ilu;

133:   ilu                = (PC_ILU*)pc->data;
134:   ilu->reuseordering = flag;
135:   return(0);
136: }
137: EXTERN_C_END

139: EXTERN_C_BEGIN
140: #undef __FUNCT__  
142: int PCILUDTSetReuseFill_ILUDT(PC pc,PetscTruth flag)
143: {
144:   PC_ILU *ilu;

147:   ilu = (PC_ILU*)pc->data;
148:   ilu->reusefill = flag;
149:   if (flag) SETERRQ(1,"Not yet supported");
150:   return(0);
151: }
152: EXTERN_C_END

154: EXTERN_C_BEGIN
155: #undef __FUNCT__  
157: int PCILUSetLevels_ILU(PC pc,int levels)
158: {
159:   PC_ILU *ilu;
160:   int    ierr;

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

165:   if (!pc->setupcalled) {
166:     ilu->info.levels = levels;
167:   } else if (ilu->usedt == PETSC_TRUE || ilu->info.levels != levels) {
168:     ilu->info.levels = levels;
169:     pc->setupcalled  = 0;
170:     ilu->usedt       = PETSC_FALSE;
171:     PCDestroy_ILU_Internal(pc);
172:   }
173:   return(0);
174: }
175: EXTERN_C_END

177: EXTERN_C_BEGIN
178: #undef __FUNCT__  
180: int PCILUSetUseInPlace_ILU(PC pc)
181: {
182:   PC_ILU *dir;

185:   dir          = (PC_ILU*)pc->data;
186:   dir->inplace = PETSC_TRUE;
187:   return(0);
188: }
189: EXTERN_C_END

191: EXTERN_C_BEGIN
192: #undef __FUNCT__  
194: int PCILUSetAllowDiagonalFill_ILU(PC pc)
195: {
196:   PC_ILU *dir;

199:   dir                 = (PC_ILU*)pc->data;
200:   dir->info.diagonal_fill = 1;
201:   return(0);
202: }
203: EXTERN_C_END

205: #undef __FUNCT__  
207: /*@
208:    PCILUSetZeroPivot - Sets the size at which smaller pivots are declared to be zero

210:    Collective on PC
211:    
212:    Input Parameters:
213: +  pc - the preconditioner context
214: -  zero - all pivots smaller than this will be considered zero

216:    Options Database Key:
217: .  -pc_ilu_zeropivot <zero> - Sets the zero pivot size

219:    Level: intermediate

221: .keywords: PC, set, factorization, direct, fill

223: .seealso: PCILUSetFill(), PCLUSetDamp(), PCLUSetZeroPivot()
224: @*/
225: int PCILUSetZeroPivot(PC pc,PetscReal zero)
226: {
227:   int ierr,(*f)(PC,PetscReal);

231:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetZeroPivot_C",(void (**)(void))&f);
232:   if (f) {
233:     (*f)(pc,zero);
234:   }
235:   return(0);
236: }

238: #undef __FUNCT__  
240: /*@
241:    PCILUSetDamping - adds this quantity to the diagonal of the matrix during the 
242:      ILU numerical factorization

244:    Collective on PC
245:    
246:    Input Parameters:
247: +  pc - the preconditioner context
248: -  damping - amount of damping

250:    Options Database Key:
251: .  -pc_ilu_damping <damping> - Sets damping amount or PETSC_DECIDE for the default

253:    Note: If 0.0 is given, then no damping is used. If a diagonal element is classified as a zero
254:          pivot, then the damping is doubled until this is alleviated.

256:    Level: intermediate

258: .keywords: PC, set, factorization, direct, fill

260: .seealso: PCILUSetFill(), PCLUSetDamp()
261: @*/
262: int PCILUSetDamping(PC pc,PetscReal damping)
263: {
264:   int ierr,(*f)(PC,PetscReal);

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

275: #undef __FUNCT__  
277: /*@
278:    PCILUSetUseDropTolerance - The preconditioner will use an ILU 
279:    based on a drop tolerance.

281:    Collective on PC

283:    Input Parameters:
284: +  pc - the preconditioner context
285: .  dt - the drop tolerance, try from 1.e-10 to .1
286: .  dtcol - tolerance for column pivot, good values [0.1 to 0.01]
287: -  maxrowcount - the max number of nonzeros allowed in a row, best value
288:                  depends on the number of nonzeros in row of original matrix

290:    Options Database Key:
291: .  -pc_ilu_use_drop_tolerance <dt,dtcol,maxrowcount> - Sets drop tolerance

293:    Level: intermediate

295:     Notes:
296:       This uses the iludt() code of Saad's SPARSKIT package

298: .keywords: PC, levels, reordering, factorization, incomplete, ILU
299: @*/
300: int PCILUSetUseDropTolerance(PC pc,PetscReal dt,PetscReal dtcol,int maxrowcount)
301: {
302:   int ierr,(*f)(PC,PetscReal,PetscReal,int);

306:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetUseDropTolerance_C",(void (**)(void))&f);
307:   if (f) {
308:     (*f)(pc,dt,dtcol,maxrowcount);
309:   }
310:   return(0);
311: }

313: #undef __FUNCT__  
315: /*@
316:    PCILUSetFill - Indicate the amount of fill you expect in the factored matrix,
317:    where fill = number nonzeros in factor/number nonzeros in original matrix.

319:    Collective on PC

321:    Input Parameters:
322: +  pc - the preconditioner context
323: -  fill - amount of expected fill

325:    Options Database Key:
326: $  -pc_ilu_fill <fill>

328:    Note:
329:    For sparse matrix factorizations it is difficult to predict how much 
330:    fill to expect. By running with the option -log_info PETSc will print the 
331:    actual amount of fill used; allowing you to set the value accurately for
332:    future runs. But default PETSc uses a value of 1.0

334:    Level: intermediate

336: .keywords: PC, set, factorization, direct, fill

338: .seealso: PCLUSetFill()
339: @*/
340: int PCILUSetFill(PC pc,PetscReal fill)
341: {
342:   int ierr,(*f)(PC,PetscReal);

346:   if (fill < 1.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Fill factor cannot be less than 1.0");
347:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetFill_C",(void (**)(void))&f);
348:   if (f) {
349:     (*f)(pc,fill);
350:   }
351:   return(0);
352: }

354: #undef __FUNCT__  
356: /*@C
357:     PCILUSetMatOrdering - Sets the ordering routine (to reduce fill) to 
358:     be used in the ILU factorization.

360:     Collective on PC

362:     Input Parameters:
363: +   pc - the preconditioner context
364: -   ordering - the matrix ordering name, for example, MATORDERING_ND or MATORDERING_RCM

366:     Options Database Key:
367: .   -pc_ilu_mat_ordering_type <nd,rcm,...> - Sets ordering routine

369:     Level: intermediate

371:     Notes: natural ordering is used by default

373: .seealso: PCLUSetMatOrdering()

375: .keywords: PC, ILU, set, matrix, reordering

377: @*/
378: int PCILUSetMatOrdering(PC pc,MatOrderingType ordering)
379: {
380:   int ierr,(*f)(PC,MatOrderingType);

384:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetMatOrdering_C",(void (**)(void))&f);
385:   if (f) {
386:     (*f)(pc,ordering);
387:   }
388:   return(0);
389: }

391: #undef __FUNCT__  
393: /*@
394:    PCILUSetReuseOrdering - When similar matrices are factored, this
395:    causes the ordering computed in the first factor to be used for all
396:    following factors; applies to both fill and drop tolerance ILUs.

398:    Collective on PC

400:    Input Parameters:
401: +  pc - the preconditioner context
402: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

404:    Options Database Key:
405: .  -pc_ilu_reuse_ordering - Activate PCILUSetReuseOrdering()

407:    Level: intermediate

409: .keywords: PC, levels, reordering, factorization, incomplete, ILU

411: .seealso: PCILUDTSetReuseFill(), PCLUSetReuseOrdering(), PCLUSetReuseFill()
412: @*/
413: int PCILUSetReuseOrdering(PC pc,PetscTruth flag)
414: {
415:   int ierr,(*f)(PC,PetscTruth);

419:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetReuseOrdering_C",(void (**)(void))&f);
420:   if (f) {
421:     (*f)(pc,flag);
422:   }
423:   return(0);
424: }

426: #undef __FUNCT__  
428: /*@
429:    PCILUDTSetReuseFill - When matrices with same nonzero structure are ILUDT factored,
430:    this causes later ones to use the fill computed in the initial factorization.

432:    Collective on PC

434:    Input Parameters:
435: +  pc - the preconditioner context
436: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

438:    Options Database Key:
439: .  -pc_iludt_reuse_fill - Activates PCILUDTSetReuseFill()

441:    Level: intermediate

443: .keywords: PC, levels, reordering, factorization, incomplete, ILU

445: .seealso: PCILUSetReuseOrdering(), PCLUSetReuseOrdering(), PCLUSetReuseFill()
446: @*/
447: int PCILUDTSetReuseFill(PC pc,PetscTruth flag)
448: {
449:   int ierr,(*f)(PC,PetscTruth);

453:   PetscObjectQueryFunction((PetscObject)pc,"PCILUDTSetReuseFill_C",(void (**)(void))&f);
454:   if (f) {
455:     (*f)(pc,flag);
456:   }
457:   return(0);
458: }

460: #undef __FUNCT__  
462: /*@
463:    PCILUSetLevels - Sets the number of levels of fill to use.

465:    Collective on PC

467:    Input Parameters:
468: +  pc - the preconditioner context
469: -  levels - number of levels of fill

471:    Options Database Key:
472: .  -pc_ilu_levels <levels> - Sets fill level

474:    Level: intermediate

476: .keywords: PC, levels, fill, factorization, incomplete, ILU
477: @*/
478: int PCILUSetLevels(PC pc,int levels)
479: {
480:   int ierr,(*f)(PC,int);

484:   if (levels < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"negative levels");
485:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetLevels_C",(void (**)(void))&f);
486:   if (f) {
487:     (*f)(pc,levels);
488:   }
489:   return(0);
490: }

492: #undef __FUNCT__  
494: /*@
495:    PCILUSetAllowDiagonalFill - Causes all diagonal matrix entries to be 
496:    treated as level 0 fill even if there is no non-zero location.

498:    Collective on PC

500:    Input Parameters:
501: +  pc - the preconditioner context

503:    Options Database Key:
504: .  -pc_ilu_diagonal_fill

506:    Notes:
507:    Does not apply with 0 fill.

509:    Level: intermediate

511: .keywords: PC, levels, fill, factorization, incomplete, ILU
512: @*/
513: int PCILUSetAllowDiagonalFill(PC pc)
514: {
515:   int ierr,(*f)(PC);

519:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetAllowDiagonalFill_C",(void (**)(void))&f);
520:   if (f) {
521:     (*f)(pc);
522:   }
523:   return(0);
524: }

526: #undef __FUNCT__  
528: /*@
529:    PCILUSetUseInPlace - Tells the system to do an in-place incomplete factorization.
530:    Collective on PC

532:    Input Parameters:
533: .  pc - the preconditioner context

535:    Options Database Key:
536: .  -pc_ilu_in_place - Activates in-place factorization

538:    Notes:
539:    PCILUSetUseInPlace() is intended for use with matrix-free variants of
540:    Krylov methods, or when a different matrices are employed for the linear
541:    system and preconditioner, or with ASM preconditioning.  Do NOT use 
542:    this option if the linear system
543:    matrix also serves as the preconditioning matrix, since the factored
544:    matrix would then overwrite the original matrix. 

546:    Only works well with ILU(0).

548:    Level: intermediate

550: .keywords: PC, set, factorization, inplace, in-place, ILU

552: .seealso:  PCLUSetUseInPlace()
553: @*/
554: int PCILUSetUseInPlace(PC pc)
555: {
556:   int ierr,(*f)(PC);

560:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetUseInPlace_C",(void (**)(void))&f);
561:   if (f) {
562:     (*f)(pc);
563:   }
564:   return(0);
565: }

567: #undef __FUNCT__  
569: /*@
570:     PCILUSetPivotInBlocks - Determines if pivoting is done while factoring each block
571:       with BAIJ or SBAIJ matrices

573:     Collective on PC

575:     Input Parameters:
576: +   pc - the preconditioner context
577: -   pivot - PETSC_TRUE or PETSC_FALSE

579:     Options Database Key:
580: .   -pc_ilu_pivot_in_blocks <true,false>

582:     Level: intermediate

584: .seealso: PCIILUSetMatOrdering(), PCILUSetPivoting()
585: @*/
586: int PCILUSetPivotInBlocks(PC pc,PetscTruth pivot)
587: {
588:   int ierr,(*f)(PC,PetscTruth);

591:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetPivotInBlocks_C",(void (**)(void))&f);
592:   if (f) {
593:     (*f)(pc,pivot);
594:   }
595:   return(0);
596: }

598: /* ------------------------------------------------------------------------------------------*/

600: EXTERN_C_BEGIN
601: #undef __FUNCT__  
603: int PCILUSetPivotInBlocks_ILU(PC pc,PetscTruth pivot)
604: {
605:   PC_ILU *dir = (PC_ILU*)pc->data;

608:   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
609:   return(0);
610: }
611: EXTERN_C_END

613: #undef __FUNCT__  
615: static int PCSetFromOptions_ILU(PC pc)
616: {
617:   int        ierr,dtmax = 3,itmp;
618:   PetscTruth flg,set;
619:   PetscReal  dt[3];
620:   char       tname[256];
621:   PC_ILU     *ilu = (PC_ILU*)pc->data;
622:   PetscFList ordlist;
623:   PetscReal  tol;

626:   MatOrderingRegisterAll(PETSC_NULL);
627:   PetscOptionsHead("ILU Options");
628:     PetscOptionsInt("-pc_ilu_levels","levels of fill","PCILUSetLevels",(int)ilu->info.levels,&itmp,&flg);
629:     if (flg) ilu->info.levels = (double) itmp;
630:     PetscOptionsName("-pc_ilu_in_place","do factorization in place","PCILUSetUseInPlace",&ilu->inplace);
631:     PetscOptionsName("-pc_ilu_diagonal_fill","Allow fill into empty diagonal entry","PCILUSetAllowDiagonalFill",&flg);
632:     ilu->info.diagonal_fill = (double) flg;
633:     PetscOptionsName("-pc_iludt_reuse_fill","Reuse fill from previous ILUdt","PCILUDTSetReuseFill",&ilu->reusefill);
634:     PetscOptionsName("-pc_ilu_reuse_ordering","Reuse previous reordering","PCILUSetReuseOrdering",&ilu->reuseordering);
635:     PetscOptionsName("-pc_ilu_damping","Damping added to diagonal","PCILUSetDamping",&flg);
636:     if (flg) {
637:       PCILUSetDamping(pc,(PetscReal) PETSC_DECIDE);
638:     }
639:     PetscOptionsReal("-pc_ilu_damping","Damping added to diagonal","PCILUSetDamping",ilu->info.damping,&ilu->info.damping,0);
640:     PetscOptionsReal("-pc_ilu_zeropivot","Pivot is considered zero if less than","PCILUSetSetZeroPivot",ilu->info.zeropivot,&ilu->info.zeropivot,0);

642:     dt[0] = ilu->info.dt;
643:     dt[1] = ilu->info.dtcol;
644:     dt[2] = ilu->info.dtcount;
645:     PetscOptionsRealArray("-pc_ilu_use_drop_tolerance","<dt,dtcol,maxrowcount>","PCILUSetUseDropTolerance",dt,&dtmax,&flg);
646:     if (flg) {
647:       PCILUSetUseDropTolerance(pc,dt[0],dt[1],(int)dt[2]);
648:     }
649:     PetscOptionsReal("-pc_ilu_fill","Expected fill in factorization","PCILUSetFill",ilu->info.fill,&ilu->info.fill,&flg);
650:     PetscOptionsReal("-pc_ilu_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","MatReorderForNonzeroDiagonal",0.0,&tol,0);

652:     MatGetOrderingList(&ordlist);
653:     PetscOptionsList("-pc_ilu_mat_ordering_type","Reorder to reduce nonzeros in ILU","PCILUSetMatOrdering",ordlist,ilu->ordering,tname,256,&flg);
654:     if (flg) {
655:       PCILUSetMatOrdering(pc,tname);
656:     }
657:     flg = ilu->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
658:     PetscOptionsLogical("-pc_ilu_pivot_in_blocks","Pivot inside matrix blocks for BAIJ and SBAIJ","PCILUSetPivotInBlocks",flg,&flg,&set);
659:     if (set) {
660:       PCILUSetPivotInBlocks(pc,flg);
661:     }
662:   PetscOptionsTail();
663:   return(0);
664: }

666: #undef __FUNCT__  
668: static int PCView_ILU(PC pc,PetscViewer viewer)
669: {
670:   PC_ILU     *ilu = (PC_ILU*)pc->data;
671:   int        ierr;
672:   PetscTruth isstring,isascii;

675:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
676:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
677:   if (isascii) {
678:     if (ilu->usedt) {
679:         PetscViewerASCIIPrintf(viewer,"  ILU: drop tolerance %gn",ilu->info.dt);
680:         PetscViewerASCIIPrintf(viewer,"  ILU: max nonzeros per row %dn",(int)ilu->info.dtcount);
681:         PetscViewerASCIIPrintf(viewer,"  ILU: column permutation tolerance %gn",ilu->info.dtcol);
682:     } else if (ilu->info.levels == 1) {
683:         PetscViewerASCIIPrintf(viewer,"  ILU: %d level of filln",(int)ilu->info.levels);
684:     } else {
685:         PetscViewerASCIIPrintf(viewer,"  ILU: %d levels of filln",(int)ilu->info.levels);
686:     }
687:     PetscViewerASCIIPrintf(viewer,"  ILU: max fill ratio allocated %gn",ilu->info.fill);
688:     PetscViewerASCIIPrintf(viewer,"  ILU: tolerance for zero pivot %gn",ilu->info.zeropivot);
689:     if (ilu->inplace) {PetscViewerASCIIPrintf(viewer,"       in-place factorizationn");}
690:     else              {PetscViewerASCIIPrintf(viewer,"       out-of-place factorizationn");}
691:     PetscViewerASCIIPrintf(viewer,"       matrix ordering: %sn",ilu->ordering);
692:     if (ilu->reusefill)     {PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorizationn");}
693:     if (ilu->reuseordering) {PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorizationn");}
694:     if (ilu->fact) {
695:       PetscViewerASCIIPrintf(viewer,"       Factored matrix followsn");
696:       PetscViewerASCIIPushTab(viewer);
697:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
698:       MatView(ilu->fact,viewer);
699:       PetscViewerPopFormat(viewer);
700:       PetscViewerASCIIPopTab(viewer);
701:     }
702:   } else if (isstring) {
703:     PetscViewerStringSPrintf(viewer," lvls=%g,order=%s",ilu->info.levels,ilu->ordering);
704:   } else {
705:     SETERRQ1(1,"Viewer type %s not supported for PCILU",((PetscObject)viewer)->type_name);
706:   }
707:   return(0);
708: }

710: #undef __FUNCT__  
712: static int PCSetUp_ILU(PC pc)
713: {
714:   int        ierr;
715:   PetscTruth flg;
716:   PC_ILU     *ilu = (PC_ILU*)pc->data;

719:   if (ilu->inplace) {
720:     if (!pc->setupcalled) {

722:       /* In-place factorization only makes sense with the natural ordering,
723:          so we only need to get the ordering once, even if nonzero structure changes */
724:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
725:       if (ilu->row) PetscLogObjectParent(pc,ilu->row);
726:       if (ilu->col) PetscLogObjectParent(pc,ilu->col);
727:     }

729:     /* In place ILU only makes sense with fill factor of 1.0 because 
730:        cannot have levels of fill */
731:     ilu->info.fill          = 1.0;
732:     ilu->info.diagonal_fill = 0;
733:     MatILUFactor(pc->pmat,ilu->row,ilu->col,&ilu->info);
734:     ilu->fact = pc->pmat;
735:   } else if (ilu->usedt) {
736:     if (!pc->setupcalled) {
737:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
738:       if (ilu->row) PetscLogObjectParent(pc,ilu->row);
739:       if (ilu->col) PetscLogObjectParent(pc,ilu->col);
740:       MatILUDTFactor(pc->pmat,&ilu->info,ilu->row,ilu->col,&ilu->fact);
741:       PetscLogObjectParent(pc,ilu->fact);
742:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
743:       MatDestroy(ilu->fact);
744:       if (!ilu->reuseordering) {
745:         if (ilu->row) {ISDestroy(ilu->row);}
746:         if (ilu->col) {ISDestroy(ilu->col);}
747:         MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
748:         if (ilu->row) PetscLogObjectParent(pc,ilu->row);
749:         if (ilu->col) PetscLogObjectParent(pc,ilu->col);
750:       }
751:       MatILUDTFactor(pc->pmat,&ilu->info,ilu->row,ilu->col,&ilu->fact);
752:       PetscLogObjectParent(pc,ilu->fact);
753:     } else if (!ilu->reusefill) {
754:       MatDestroy(ilu->fact);
755:       MatILUDTFactor(pc->pmat,&ilu->info,ilu->row,ilu->col,&ilu->fact);
756:       PetscLogObjectParent(pc,ilu->fact);
757:     } else {
758:       MatLUFactorNumeric(pc->pmat,&ilu->fact);
759:     }
760:    } else {
761:     if (!pc->setupcalled) {
762:       /* first time in so compute reordering and symbolic factorization */
763:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
764:       if (ilu->row) PetscLogObjectParent(pc,ilu->row);
765:       if (ilu->col) PetscLogObjectParent(pc,ilu->col);
766:       /*  Remove zeros along diagonal?     */
767:       PetscOptionsHasName(pc->prefix,"-pc_ilu_nonzeros_along_diagonal",&flg);
768:       if (flg) {
769:         PetscReal ntol = 1.e-10;
770:         PetscOptionsGetReal(pc->prefix,"-pc_ilu_nonzeros_along_diagonal",&ntol,PETSC_NULL);
771:         MatReorderForNonzeroDiagonal(pc->pmat,ntol,ilu->row,ilu->col);
772:       }

774:       MatILUFactorSymbolic(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
775:       PetscLogObjectParent(pc,ilu->fact);
776:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
777:       if (!ilu->reuseordering) {
778:         /* compute a new ordering for the ILU */
779:         ISDestroy(ilu->row);
780:         ISDestroy(ilu->col);
781:         MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
782:         if (ilu->row) PetscLogObjectParent(pc,ilu->row);
783:         if (ilu->col) PetscLogObjectParent(pc,ilu->col);
784:         /*  Remove zeros along diagonal?     */
785:         PetscOptionsHasName(pc->prefix,"-pc_ilu_nonzeros_along_diagonal",&flg);
786:         if (flg) {
787:           PetscReal ntol = 1.e-10;
788:           PetscOptionsGetReal(pc->prefix,"-pc_ilu_nonzeros_along_diagonal",&ntol,PETSC_NULL);
789:           MatReorderForNonzeroDiagonal(pc->pmat,ntol,ilu->row,ilu->col);
790:         }
791:       }
792:       MatDestroy(ilu->fact);
793:       MatILUFactorSymbolic(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
794:       PetscLogObjectParent(pc,ilu->fact);
795:     }
796:     MatLUFactorNumeric(pc->pmat,&ilu->fact);
797:   }
798:   return(0);
799: }

801: #undef __FUNCT__  
803: static int PCDestroy_ILU(PC pc)
804: {
805:   PC_ILU *ilu = (PC_ILU*)pc->data;
806:   int    ierr;

809:   PCDestroy_ILU_Internal(pc);
810:   PetscStrfree(ilu->ordering);
811:   PetscFree(ilu);
812:   return(0);
813: }

815: #undef __FUNCT__  
817: static int PCApply_ILU(PC pc,Vec x,Vec y)
818: {
819:   PC_ILU *ilu = (PC_ILU*)pc->data;
820:   int    ierr;

823:   MatSolve(ilu->fact,x,y);
824:   return(0);
825: }

827: #undef __FUNCT__  
829: static int PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
830: {
831:   PC_ILU *ilu = (PC_ILU*)pc->data;
832:   int    ierr;

835:   MatSolveTranspose(ilu->fact,x,y);
836:   return(0);
837: }

839: #undef __FUNCT__  
841: static int PCGetFactoredMatrix_ILU(PC pc,Mat *mat)
842: {
843:   PC_ILU *ilu = (PC_ILU*)pc->data;

846:   if (!ilu->fact) SETERRQ(1,"Matrix not yet factored; call after SLESSetUp() or PCSetUp()");
847:   *mat = ilu->fact;
848:   return(0);
849: }

851: EXTERN_C_BEGIN
852: #undef __FUNCT__  
854: int PCCreate_ILU(PC pc)
855: {
856:   int    ierr;
857:   PC_ILU *ilu;

860:   PetscNew(PC_ILU,&ilu);
861:   PetscLogObjectMemory(pc,sizeof(PC_ILU));

863:   ilu->fact                    = 0;
864:   ilu->info.levels             = 0;
865:   ilu->info.fill               = 1.0;
866:   ilu->col                     = 0;
867:   ilu->row                     = 0;
868:   ilu->inplace                 = PETSC_FALSE;
869:   PetscStrallocpy(MATORDERING_NATURAL,&ilu->ordering);
870:   ilu->reuseordering           = PETSC_FALSE;
871:   ilu->usedt                   = PETSC_FALSE;
872:   ilu->info.dt                 = PETSC_DEFAULT;
873:   ilu->info.dtcount            = PETSC_DEFAULT;
874:   ilu->info.dtcol              = PETSC_DEFAULT;
875:   ilu->info.damping            = 0.0;
876:   ilu->info.zeropivot          = 1.e-12;
877:   ilu->info.pivotinblocks      = 1.0;
878:   ilu->reusefill               = PETSC_FALSE;
879:   ilu->info.diagonal_fill      = 0;
880:   pc->data                     = (void*)ilu;

882:   pc->ops->destroy             = PCDestroy_ILU;
883:   pc->ops->apply               = PCApply_ILU;
884:   pc->ops->applytranspose      = PCApplyTranspose_ILU;
885:   pc->ops->setup               = PCSetUp_ILU;
886:   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
887:   pc->ops->getfactoredmatrix   = PCGetFactoredMatrix_ILU;
888:   pc->ops->view                = PCView_ILU;
889:   pc->ops->applyrichardson     = 0;

891:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetUseDropTolerance_C","PCILUSetUseDropTolerance_ILU",
892:                     PCILUSetUseDropTolerance_ILU);
893:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetFill_C","PCILUSetFill_ILU",
894:                     PCILUSetFill_ILU);
895:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetDamping_C","PCILUSetDamping_ILU",
896:                     PCILUSetDamping_ILU);
897:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetMatOrdering_C","PCILUSetMatOrdering_ILU",
898:                     PCILUSetMatOrdering_ILU);
899:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetReuseOrdering_C","PCILUSetReuseOrdering_ILU",
900:                     PCILUSetReuseOrdering_ILU);
901:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUDTSetReuseFill_C","PCILUDTSetReuseFill_ILUDT",
902:                     PCILUDTSetReuseFill_ILUDT);
903:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetLevels_C","PCILUSetLevels_ILU",
904:                     PCILUSetLevels_ILU);
905:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetUseInPlace_C","PCILUSetUseInPlace_ILU",
906:                     PCILUSetUseInPlace_ILU);
907:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetAllowDiagonalFill_C","PCILUSetAllowDiagonalFill_ILU",
908:                     PCILUSetAllowDiagonalFill_ILU);
909:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetPivotInBlocks_C","PCILUSetPivotInBlocks_ILU",
910:                     PCILUSetPivotInBlocks_ILU);
911:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetZeroPivot_C","PCILUSetZeroPivot_ILU",
912:                     PCILUSetZeroPivot_ILU);
913:   return(0);
914: }
915: EXTERN_C_END