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