xref: /petsc/src/ksp/pc/impls/factor/ilu/ilu.c (revision ce94432eddcd14845bc7e8083b7f8ea723b9bf7d)
1dba47a55SKris Buschelman 
29b54502bSHong Zhang /*
39b54502bSHong Zhang    Defines a ILU factorization preconditioner for any Mat implementation
49b54502bSHong Zhang */
5c6db04a5SJed Brown #include <../src/ksp/pc/impls/factor/ilu/ilu.h>     /*I "petscpc.h"  I*/
69b54502bSHong Zhang 
79b54502bSHong Zhang /* ------------------------------------------------------------------------------------------*/
8afaefe49SHong Zhang EXTERN_C_BEGIN
9afaefe49SHong Zhang #undef __FUNCT__
102401956bSBarry Smith #define __FUNCT__ "PCFactorSetReuseFill_ILU"
117087cfbeSBarry Smith PetscErrorCode  PCFactorSetReuseFill_ILU(PC pc,PetscBool flag)
122401956bSBarry Smith {
13075768bcSBarry Smith   PC_ILU *lu = (PC_ILU*)pc->data;
142401956bSBarry Smith 
152401956bSBarry Smith   PetscFunctionBegin;
162401956bSBarry Smith   lu->reusefill = flag;
172401956bSBarry Smith   PetscFunctionReturn(0);
182401956bSBarry Smith }
192401956bSBarry Smith EXTERN_C_END
202401956bSBarry Smith 
212401956bSBarry Smith EXTERN_C_BEGIN
222401956bSBarry Smith #undef __FUNCT__
232401956bSBarry Smith #define __FUNCT__ "PCFactorReorderForNonzeroDiagonal_ILU"
247087cfbeSBarry Smith PetscErrorCode  PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)
259b54502bSHong Zhang {
269b54502bSHong Zhang   PC_ILU *ilu = (PC_ILU*)pc->data;
279b54502bSHong Zhang 
289b54502bSHong Zhang   PetscFunctionBegin;
299b54502bSHong Zhang   ilu->nonzerosalongdiagonal = PETSC_TRUE;
302fa5cd67SKarl Rupp   if (z == PETSC_DECIDE) ilu->nonzerosalongdiagonaltol = 1.e-10;
312fa5cd67SKarl Rupp   else ilu->nonzerosalongdiagonaltol = z;
329b54502bSHong Zhang   PetscFunctionReturn(0);
339b54502bSHong Zhang }
349b54502bSHong Zhang EXTERN_C_END
359b54502bSHong Zhang 
369b54502bSHong Zhang #undef __FUNCT__
37574deadeSBarry Smith #define __FUNCT__ "PCReset_ILU"
38574deadeSBarry Smith PetscErrorCode PCReset_ILU(PC pc)
399b54502bSHong Zhang {
409b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
419b54502bSHong Zhang   PetscErrorCode ierr;
429b54502bSHong Zhang 
439b54502bSHong Zhang   PetscFunctionBegin;
44298cc208SBarry Smith   if (!ilu->inplace) {ierr = MatDestroy(&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);}
456bf464f9SBarry Smith   if (ilu->row && ilu->col && ilu->row != ilu->col) {ierr = ISDestroy(&ilu->row);CHKERRQ(ierr);}
46fcfd50ebSBarry Smith   ierr = ISDestroy(&ilu->col);CHKERRQ(ierr);
479b54502bSHong Zhang   PetscFunctionReturn(0);
489b54502bSHong Zhang }
499b54502bSHong Zhang 
509b54502bSHong Zhang EXTERN_C_BEGIN
519b54502bSHong Zhang #undef __FUNCT__
52b7c853c4SBarry Smith #define __FUNCT__ "PCFactorSetDropTolerance_ILU"
537087cfbeSBarry Smith PetscErrorCode  PCFactorSetDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
549b54502bSHong Zhang {
55075768bcSBarry Smith   PC_ILU *ilu = (PC_ILU*)pc->data;
569b54502bSHong Zhang 
579b54502bSHong Zhang   PetscFunctionBegin;
584c9036c7SBarry Smith   if (pc->setupcalled && (((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
59*ce94432eSBarry Smith     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot change drop tolerance after using PC");
609b54502bSHong Zhang   }
61075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dt      = dt;
62075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcol   = dtcol;
63075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcount = dtcount;
644c9036c7SBarry Smith   ((PC_Factor*)ilu)->info.usedt   = 1.0;
659b54502bSHong Zhang   PetscFunctionReturn(0);
669b54502bSHong Zhang }
679b54502bSHong Zhang EXTERN_C_END
689b54502bSHong Zhang 
699b54502bSHong Zhang EXTERN_C_BEGIN
709b54502bSHong Zhang #undef __FUNCT__
712401956bSBarry Smith #define __FUNCT__ "PCFactorSetReuseOrdering_ILU"
727087cfbeSBarry Smith PetscErrorCode  PCFactorSetReuseOrdering_ILU(PC pc,PetscBool flag)
739b54502bSHong Zhang {
74075768bcSBarry Smith   PC_ILU *ilu = (PC_ILU*)pc->data;
759b54502bSHong Zhang 
769b54502bSHong Zhang   PetscFunctionBegin;
779b54502bSHong Zhang   ilu->reuseordering = flag;
789b54502bSHong Zhang   PetscFunctionReturn(0);
799b54502bSHong Zhang }
809b54502bSHong Zhang EXTERN_C_END
819b54502bSHong Zhang 
829b54502bSHong Zhang EXTERN_C_BEGIN
839b54502bSHong Zhang #undef __FUNCT__
842401956bSBarry Smith #define __FUNCT__ "PCFactorSetUseInPlace_ILU"
857087cfbeSBarry Smith PetscErrorCode  PCFactorSetUseInPlace_ILU(PC pc)
869b54502bSHong Zhang {
87075768bcSBarry Smith   PC_ILU *dir = (PC_ILU*)pc->data;
889b54502bSHong Zhang 
899b54502bSHong Zhang   PetscFunctionBegin;
909b54502bSHong Zhang   dir->inplace = PETSC_TRUE;
919b54502bSHong Zhang   PetscFunctionReturn(0);
929b54502bSHong Zhang }
939b54502bSHong Zhang EXTERN_C_END
949b54502bSHong Zhang 
959b54502bSHong Zhang #undef __FUNCT__
969b54502bSHong Zhang #define __FUNCT__ "PCSetFromOptions_ILU"
979b54502bSHong Zhang static PetscErrorCode PCSetFromOptions_ILU(PC pc)
989b54502bSHong Zhang {
999b54502bSHong Zhang   PetscErrorCode ierr;
10078fc6b22SHong Zhang   PetscInt       itmp;
101ace3abfcSBarry Smith   PetscBool      flg;
1029b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
1039b54502bSHong Zhang   PetscReal      tol;
1042fa5cd67SKarl Rupp   /* PetscReal      dt[3]; */
1059b54502bSHong Zhang 
1069b54502bSHong Zhang   PetscFunctionBegin;
1079b54502bSHong Zhang   ierr = PetscOptionsHead("ILU Options");CHKERRQ(ierr);
1088ff23777SHong Zhang   ierr = PCSetFromOptions_Factor(pc);CHKERRQ(ierr);
1098ff23777SHong Zhang 
110075768bcSBarry Smith   ierr = PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);CHKERRQ(ierr);
111075768bcSBarry Smith   if (flg) ((PC_Factor*)ilu)->info.levels = itmp;
1122fa5cd67SKarl Rupp 
11390d69ab7SBarry Smith   flg  = PETSC_FALSE;
1140298fd71SBarry Smith   ierr = PetscOptionsBool("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",flg,&flg,NULL);CHKERRQ(ierr);
115075768bcSBarry Smith   ((PC_Factor*)ilu)->info.diagonal_fill = (double) flg;
1165a586d82SBarry Smith   /*
117075768bcSBarry Smith   dt[0] = ((PC_Factor*)ilu)->info.dt;
118075768bcSBarry Smith   dt[1] = ((PC_Factor*)ilu)->info.dtcol;
119075768bcSBarry Smith   dt[2] = ((PC_Factor*)ilu)->info.dtcount;
1205a586d82SBarry Smith 
12178fc6b22SHong Zhang   PetscInt       dtmax = 3;
12278fc6b22SHong Zhang   ierr = PetscOptionsRealArray("-pc_factor_drop_tolerance,","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);CHKERRQ(ierr);
1239b54502bSHong Zhang   if (flg) {
124b7c853c4SBarry Smith     ierr = PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);CHKERRQ(ierr);
1259b54502bSHong Zhang   }
12678fc6b22SHong Zhang   */
1272401956bSBarry Smith   ierr = PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);CHKERRQ(ierr);
1289b54502bSHong Zhang   if (flg) {
1299b54502bSHong Zhang     tol  = PETSC_DECIDE;
1302401956bSBarry Smith     ierr = PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);CHKERRQ(ierr);
1312401956bSBarry Smith     ierr = PCFactorReorderForNonzeroDiagonal(pc,tol);CHKERRQ(ierr);
1329b54502bSHong Zhang   }
1339b54502bSHong Zhang 
1349b54502bSHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
1359b54502bSHong Zhang   PetscFunctionReturn(0);
1369b54502bSHong Zhang }
1379b54502bSHong Zhang 
1389b54502bSHong Zhang #undef __FUNCT__
1399b54502bSHong Zhang #define __FUNCT__ "PCView_ILU"
1409b54502bSHong Zhang static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
1419b54502bSHong Zhang {
1429b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
1439b54502bSHong Zhang   PetscErrorCode ierr;
144ace3abfcSBarry Smith   PetscBool      iascii;
1459b54502bSHong Zhang 
1469b54502bSHong Zhang   PetscFunctionBegin;
147251f4c67SDmitry Karpeev   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1489b54502bSHong Zhang   if (iascii) {
149914a5d51SHong Zhang     if (ilu->inplace) {
150914a5d51SHong Zhang       ierr = PetscViewerASCIIPrintf(viewer,"  ILU: in-place factorization\n");CHKERRQ(ierr);
1519b54502bSHong Zhang     } else {
152914a5d51SHong Zhang       ierr = PetscViewerASCIIPrintf(viewer,"  ILU: out-of-place factorization\n");CHKERRQ(ierr);
153145b9266SHong Zhang     }
154145b9266SHong Zhang 
155914a5d51SHong Zhang     if (ilu->reusefill)     {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: Reusing fill from past factorization\n");CHKERRQ(ierr);}
156914a5d51SHong Zhang     if (ilu->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: Reusing reordering from past factorization\n");CHKERRQ(ierr);}
1579b54502bSHong Zhang   }
158914a5d51SHong Zhang   ierr = PCView_Factor(pc,viewer);CHKERRQ(ierr);
1599b54502bSHong Zhang   PetscFunctionReturn(0);
1609b54502bSHong Zhang }
1619b54502bSHong Zhang 
1629b54502bSHong Zhang #undef __FUNCT__
1639b54502bSHong Zhang #define __FUNCT__ "PCSetUp_ILU"
1649b54502bSHong Zhang static PetscErrorCode PCSetUp_ILU(PC pc)
1659b54502bSHong Zhang {
1669b54502bSHong Zhang   PetscErrorCode ierr;
1679b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
168f3a39becSBarry Smith   MatInfo        info;
169ace3abfcSBarry Smith   PetscBool      flg;
1709b54502bSHong Zhang 
1719b54502bSHong Zhang   PetscFunctionBegin;
17292927226SBarry Smith   /* ugly hack to change default, since it is not support by some matrix types */
17392927226SBarry Smith   if (((PC_Factor*)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
17422d28d08SBarry Smith     ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);CHKERRQ(ierr);
17592927226SBarry Smith     if (!flg) {
17622d28d08SBarry Smith       ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&flg);CHKERRQ(ierr);
17792927226SBarry Smith       if (!flg) {
17892927226SBarry Smith         ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
17922d28d08SBarry Smith         PetscInfo(pc,"Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO");CHKERRQ(ierr);
18092927226SBarry Smith       }
18192927226SBarry Smith     }
18292927226SBarry Smith   }
18392927226SBarry Smith 
1849b54502bSHong Zhang   if (ilu->inplace) {
185075768bcSBarry Smith     CHKMEMQ;
1869b54502bSHong Zhang     if (!pc->setupcalled) {
1879b54502bSHong Zhang 
1889b54502bSHong Zhang       /* In-place factorization only makes sense with the natural ordering,
1899b54502bSHong Zhang          so we only need to get the ordering once, even if nonzero structure changes */
190075768bcSBarry Smith       ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
191efee365bSSatish Balay       if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
192efee365bSSatish Balay       if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
1939b54502bSHong Zhang     }
1949b54502bSHong Zhang 
1959b54502bSHong Zhang     /* In place ILU only makes sense with fill factor of 1.0 because
1969b54502bSHong Zhang        cannot have levels of fill */
197075768bcSBarry Smith     ((PC_Factor*)ilu)->info.fill          = 1.0;
19875567043SBarry Smith     ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
1992fa5cd67SKarl Rupp 
200fe97e370SBarry Smith     ierr = MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);CHKMEMQ;
201075768bcSBarry Smith     ((PC_Factor*)ilu)->fact = pc->pmat;
2029b54502bSHong Zhang   } else {
2039b54502bSHong Zhang     if (!pc->setupcalled) {
2049b54502bSHong Zhang       /* first time in so compute reordering and symbolic factorization */
205075768bcSBarry Smith       ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
206efee365bSSatish Balay       if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
207efee365bSSatish Balay       if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
2089b54502bSHong Zhang       /*  Remove zeros along diagonal?     */
2099b54502bSHong Zhang       if (ilu->nonzerosalongdiagonal) {
2109b54502bSHong Zhang         ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr);
2119b54502bSHong Zhang       }
21286daa42cSBarry Smith       if (!((PC_Factor*)ilu)->fact) {
21311bfe483SHong Zhang         ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
214a1f19f5aSHong Zhang       }
215075768bcSBarry Smith       ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
216075768bcSBarry Smith       ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
2172fa5cd67SKarl Rupp 
218f3a39becSBarry Smith       ilu->actualfill = info.fill_ratio_needed;
2192fa5cd67SKarl Rupp 
220075768bcSBarry Smith       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
2219b54502bSHong Zhang     } else if (pc->flag != SAME_NONZERO_PATTERN) {
2229b54502bSHong Zhang       if (!ilu->reuseordering) {
2239b54502bSHong Zhang         /* compute a new ordering for the ILU */
224fcfd50ebSBarry Smith         ierr = ISDestroy(&ilu->row);CHKERRQ(ierr);
225fcfd50ebSBarry Smith         ierr = ISDestroy(&ilu->col);CHKERRQ(ierr);
226075768bcSBarry Smith         ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
227efee365bSSatish Balay         if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
228efee365bSSatish Balay         if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
2299b54502bSHong Zhang         /*  Remove zeros along diagonal?     */
2309b54502bSHong Zhang         if (ilu->nonzerosalongdiagonal) {
2319b54502bSHong Zhang           ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr);
2329b54502bSHong Zhang         }
2339b54502bSHong Zhang       }
2346bf464f9SBarry Smith       ierr = MatDestroy(&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
23514798fb4SJed Brown       ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
236075768bcSBarry Smith       ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
237075768bcSBarry Smith       ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
2382fa5cd67SKarl Rupp 
239f3a39becSBarry Smith       ilu->actualfill = info.fill_ratio_needed;
2402fa5cd67SKarl Rupp 
241075768bcSBarry Smith       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
2429b54502bSHong Zhang     }
243075768bcSBarry Smith     CHKMEMQ;
244075768bcSBarry Smith     ierr = MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
245075768bcSBarry Smith     CHKMEMQ;
2469b54502bSHong Zhang   }
2479b54502bSHong Zhang   PetscFunctionReturn(0);
2489b54502bSHong Zhang }
2499b54502bSHong Zhang 
2509b54502bSHong Zhang #undef __FUNCT__
2519b54502bSHong Zhang #define __FUNCT__ "PCDestroy_ILU"
2529b54502bSHong Zhang static PetscErrorCode PCDestroy_ILU(PC pc)
2539b54502bSHong Zhang {
2549b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2559b54502bSHong Zhang   PetscErrorCode ierr;
2569b54502bSHong Zhang 
2579b54502bSHong Zhang   PetscFunctionBegin;
258574deadeSBarry Smith   ierr = PCReset_ILU(pc);CHKERRQ(ierr);
259503cfb0cSBarry Smith   ierr = PetscFree(((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr);
260503cfb0cSBarry Smith   ierr = PetscFree(((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
261c31cb41cSBarry Smith   ierr = PetscFree(pc->data);CHKERRQ(ierr);
2629b54502bSHong Zhang   PetscFunctionReturn(0);
2639b54502bSHong Zhang }
2649b54502bSHong Zhang 
2659b54502bSHong Zhang #undef __FUNCT__
2669b54502bSHong Zhang #define __FUNCT__ "PCApply_ILU"
2679b54502bSHong Zhang static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
2689b54502bSHong Zhang {
2699b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2709b54502bSHong Zhang   PetscErrorCode ierr;
2719b54502bSHong Zhang 
2729b54502bSHong Zhang   PetscFunctionBegin;
273075768bcSBarry Smith   ierr = MatSolve(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr);
2749b54502bSHong Zhang   PetscFunctionReturn(0);
2759b54502bSHong Zhang }
2769b54502bSHong Zhang 
2779b54502bSHong Zhang #undef __FUNCT__
2789b54502bSHong Zhang #define __FUNCT__ "PCApplyTranspose_ILU"
2799b54502bSHong Zhang static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
2809b54502bSHong Zhang {
2819b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2829b54502bSHong Zhang   PetscErrorCode ierr;
2839b54502bSHong Zhang 
2849b54502bSHong Zhang   PetscFunctionBegin;
285075768bcSBarry Smith   ierr = MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr);
2869b54502bSHong Zhang   PetscFunctionReturn(0);
2879b54502bSHong Zhang }
2889b54502bSHong Zhang 
289f0b9ad6cSBarry Smith #undef __FUNCT__
290f0b9ad6cSBarry Smith #define __FUNCT__ "PCApplySymmetricLeft_ILU"
291f0b9ad6cSBarry Smith static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
292f0b9ad6cSBarry Smith {
293f0b9ad6cSBarry Smith   PetscErrorCode ierr;
294f0b9ad6cSBarry Smith   PC_ILU         *icc = (PC_ILU*)pc->data;
295f0b9ad6cSBarry Smith 
296f0b9ad6cSBarry Smith   PetscFunctionBegin;
297f0b9ad6cSBarry Smith   ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr);
298f0b9ad6cSBarry Smith   PetscFunctionReturn(0);
299f0b9ad6cSBarry Smith }
300f0b9ad6cSBarry Smith 
301f0b9ad6cSBarry Smith #undef __FUNCT__
302f0b9ad6cSBarry Smith #define __FUNCT__ "PCApplySymmetricRight_ILU"
303f0b9ad6cSBarry Smith static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
304f0b9ad6cSBarry Smith {
305f0b9ad6cSBarry Smith   PetscErrorCode ierr;
306f0b9ad6cSBarry Smith   PC_ILU         *icc = (PC_ILU*)pc->data;
307f0b9ad6cSBarry Smith 
308f0b9ad6cSBarry Smith   PetscFunctionBegin;
309f0b9ad6cSBarry Smith   ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr);
310f0b9ad6cSBarry Smith   PetscFunctionReturn(0);
311f0b9ad6cSBarry Smith }
312f0b9ad6cSBarry Smith 
3139b54502bSHong Zhang /*MC
3149b54502bSHong Zhang      PCILU - Incomplete factorization preconditioners.
3159b54502bSHong Zhang 
3169b54502bSHong Zhang    Options Database Keys:
3172401956bSBarry Smith +  -pc_factor_levels <k> - number of levels of fill for ILU(k)
3182401956bSBarry Smith .  -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
3199b54502bSHong Zhang                       its factorization (overwrites original matrix)
3202401956bSBarry Smith .  -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
3212401956bSBarry Smith .  -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
32255ba2a51SBarry Smith .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
3232401956bSBarry Smith .  -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
3249b54502bSHong Zhang                                    this decreases the chance of getting a zero pivot
3252401956bSBarry Smith .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
3262401956bSBarry Smith .  -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
3279b54502bSHong Zhang                              than 1 the diagonal blocks are factored with partial pivoting (this increases the
3289b54502bSHong Zhang                              stability of the ILU factorization
3299b54502bSHong Zhang 
3309b54502bSHong Zhang    Level: beginner
3319b54502bSHong Zhang 
3329b54502bSHong Zhang   Concepts: incomplete factorization
3339b54502bSHong Zhang 
334d9ba8547SBarry Smith    Notes: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)
3359b54502bSHong Zhang 
3369b54502bSHong Zhang           For BAIJ matrices this implements a point block ILU
3379b54502bSHong Zhang 
338f0b9ad6cSBarry Smith           The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
339f0b9ad6cSBarry Smith           even when the matrix is not symmetric since the U stores the diagonals of the factorization.
340f0b9ad6cSBarry Smith 
341ac793be5SBarry Smith           If you are using MATSEQAIJCUSP matrices (or MATMPIAIJCUSP matrices with block Jacobi) you must have ./configured
342ac793be5SBarry Smith           PETSc with also --download-txpetscgpu to have the triangular solves performed on the GPU (factorization is never
343ac793be5SBarry Smith           done on the GPU).
344ac793be5SBarry Smith 
345c582cd25SBarry Smith    References:
346c582cd25SBarry Smith    T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
347c582cd25SBarry Smith    self-adjoint elliptic difference equations. SIAM J. Numer. Anal., 5:559--573, 1968.
348c582cd25SBarry Smith 
349c582cd25SBarry Smith    T.A. Oliphant. An implicit numerical method for solving two-dimensional time-dependent dif-
350c582cd25SBarry Smith    fusion problems. Quart. Appl. Math., 19:221--229, 1961.
351c582cd25SBarry Smith 
352c582cd25SBarry Smith    Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST
353c582cd25SBarry Smith       http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical
354c582cd25SBarry Smith       Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
355c582cd25SBarry Smith       Science and Engineering, Kluwer, pp. 167--202.
356c582cd25SBarry Smith 
357c582cd25SBarry Smith 
3589b54502bSHong Zhang .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
359145b9266SHong Zhang            PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
360b7c853c4SBarry Smith            PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
3618ff23777SHong Zhang            PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks()
3629b54502bSHong Zhang 
3639b54502bSHong Zhang M*/
3649b54502bSHong Zhang 
3659b54502bSHong Zhang EXTERN_C_BEGIN
3669b54502bSHong Zhang #undef __FUNCT__
3679b54502bSHong Zhang #define __FUNCT__ "PCCreate_ILU"
3687087cfbeSBarry Smith PetscErrorCode  PCCreate_ILU(PC pc)
3699b54502bSHong Zhang {
3709b54502bSHong Zhang   PetscErrorCode ierr;
3719b54502bSHong Zhang   PC_ILU         *ilu;
3729b54502bSHong Zhang 
3739b54502bSHong Zhang   PetscFunctionBegin;
37438f2d2fdSLisandro Dalcin   ierr = PetscNewLog(pc,PC_ILU,&ilu);CHKERRQ(ierr);
3759b54502bSHong Zhang 
376075768bcSBarry Smith   ((PC_Factor*)ilu)->fact               = 0;
377075768bcSBarry Smith   ierr                                  = MatFactorInfoInitialize(&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
378879e8a4dSBarry Smith   ((PC_Factor*)ilu)->factortype         = MAT_FACTOR_ILU;
37975567043SBarry Smith   ((PC_Factor*)ilu)->info.levels        = 0.;
380075768bcSBarry Smith   ((PC_Factor*)ilu)->info.fill          = 1.0;
3819b54502bSHong Zhang   ilu->col                              = 0;
3829b54502bSHong Zhang   ilu->row                              = 0;
3839b54502bSHong Zhang   ilu->inplace                          = PETSC_FALSE;
3842692d6eeSBarry Smith   ierr                                  = PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr);
38519fd82e9SBarry Smith   ierr                                  = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
3869b54502bSHong Zhang   ilu->reuseordering                    = PETSC_FALSE;
387075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dt            = PETSC_DEFAULT;
388075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcount       = PETSC_DEFAULT;
389075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcol         = PETSC_DEFAULT;
3909dbfc187SHong Zhang   ((PC_Factor*)ilu)->info.shifttype     = (PetscReal)MAT_SHIFT_INBLOCKS;
3910ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.shiftamount   = 100.0*PETSC_MACHINE_EPSILON;
3920ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.zeropivot     = 100.0*PETSC_MACHINE_EPSILON;
393075768bcSBarry Smith   ((PC_Factor*)ilu)->info.pivotinblocks = 1.0;
3949b54502bSHong Zhang   ilu->reusefill                        = PETSC_FALSE;
3950ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
3969b54502bSHong Zhang   pc->data                              = (void*)ilu;
3979b54502bSHong Zhang 
398574deadeSBarry Smith   pc->ops->reset               = PCReset_ILU;
3999b54502bSHong Zhang   pc->ops->destroy             = PCDestroy_ILU;
4009b54502bSHong Zhang   pc->ops->apply               = PCApply_ILU;
4019b54502bSHong Zhang   pc->ops->applytranspose      = PCApplyTranspose_ILU;
4029b54502bSHong Zhang   pc->ops->setup               = PCSetUp_ILU;
4039b54502bSHong Zhang   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
40485317021SBarry Smith   pc->ops->getfactoredmatrix   = PCFactorGetMatrix_Factor;
4059b54502bSHong Zhang   pc->ops->view                = PCView_ILU;
406f0b9ad6cSBarry Smith   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ILU;
407f0b9ad6cSBarry Smith   pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
4089b54502bSHong Zhang   pc->ops->applyrichardson     = 0;
4099b54502bSHong Zhang 
41085317021SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",
41185317021SBarry Smith                                            PCFactorSetZeroPivot_Factor);CHKERRQ(ierr);
412d90ac83dSHong Zhang   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor",
413d90ac83dSHong Zhang                                            PCFactorSetShiftType_Factor);CHKERRQ(ierr);
414d90ac83dSHong Zhang   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor",
415d90ac83dSHong Zhang                                            PCFactorSetShiftAmount_Factor);CHKERRQ(ierr);
4167112b564SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",
4177112b564SBarry Smith                                            PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr);
41811bfe483SHong Zhang   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",
41911bfe483SHong Zhang                                            PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr);
420f8260c8fSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor",
421f8260c8fSBarry Smith                                            PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr);
422b7c853c4SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetDropTolerance_C","PCFactorSetDropTolerance_ILU",
423b7c853c4SBarry Smith                                            PCFactorSetDropTolerance_ILU);CHKERRQ(ierr);
42485317021SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",
42585317021SBarry Smith                                            PCFactorSetFill_Factor);CHKERRQ(ierr);
42685317021SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",
42785317021SBarry Smith                                            PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr);
4282401956bSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_ILU",
4292401956bSBarry Smith                                            PCFactorSetReuseOrdering_ILU);CHKERRQ(ierr);
4302401956bSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_ILU",
4312401956bSBarry Smith                                            PCFactorSetReuseFill_ILU);CHKERRQ(ierr);
43285317021SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor",
43385317021SBarry Smith                                            PCFactorSetLevels_Factor);CHKERRQ(ierr);
4342401956bSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_ILU",
4352401956bSBarry Smith                                            PCFactorSetUseInPlace_ILU);CHKERRQ(ierr);
43685317021SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C","PCFactorSetAllowDiagonalFill_Factor",
43785317021SBarry Smith                                            PCFactorSetAllowDiagonalFill_Factor);CHKERRQ(ierr);
43885317021SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetPivotInBlocks_C","PCFactorSetPivotInBlocks_Factor",
43985317021SBarry Smith                                            PCFactorSetPivotInBlocks_Factor);CHKERRQ(ierr);
4402401956bSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C","PCFactorReorderForNonzeroDiagonal_ILU",
4412401956bSBarry Smith                                            PCFactorReorderForNonzeroDiagonal_ILU);CHKERRQ(ierr);
4429b54502bSHong Zhang   PetscFunctionReturn(0);
4439b54502bSHong Zhang }
4449b54502bSHong Zhang EXTERN_C_END
445