xref: /petsc/src/ksp/pc/impls/factor/ilu/ilu.c (revision 00e125f844f1e8da6dd97b92f2bb419dde0a009f)
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 #undef __FUNCT__
92401956bSBarry Smith #define __FUNCT__ "PCFactorSetReuseFill_ILU"
107087cfbeSBarry Smith PetscErrorCode  PCFactorSetReuseFill_ILU(PC pc,PetscBool flag)
112401956bSBarry Smith {
12075768bcSBarry Smith   PC_ILU *lu = (PC_ILU*)pc->data;
132401956bSBarry Smith 
142401956bSBarry Smith   PetscFunctionBegin;
152401956bSBarry Smith   lu->reusefill = flag;
162401956bSBarry Smith   PetscFunctionReturn(0);
172401956bSBarry Smith }
182401956bSBarry Smith 
192401956bSBarry Smith #undef __FUNCT__
202401956bSBarry Smith #define __FUNCT__ "PCFactorReorderForNonzeroDiagonal_ILU"
217087cfbeSBarry Smith PetscErrorCode  PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)
229b54502bSHong Zhang {
239b54502bSHong Zhang   PC_ILU *ilu = (PC_ILU*)pc->data;
249b54502bSHong Zhang 
259b54502bSHong Zhang   PetscFunctionBegin;
269b54502bSHong Zhang   ilu->nonzerosalongdiagonal = PETSC_TRUE;
272fa5cd67SKarl Rupp   if (z == PETSC_DECIDE) ilu->nonzerosalongdiagonaltol = 1.e-10;
282fa5cd67SKarl Rupp   else ilu->nonzerosalongdiagonaltol = z;
299b54502bSHong Zhang   PetscFunctionReturn(0);
309b54502bSHong Zhang }
319b54502bSHong Zhang 
329b54502bSHong Zhang #undef __FUNCT__
33574deadeSBarry Smith #define __FUNCT__ "PCReset_ILU"
34574deadeSBarry Smith PetscErrorCode PCReset_ILU(PC pc)
359b54502bSHong Zhang {
369b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
379b54502bSHong Zhang   PetscErrorCode ierr;
389b54502bSHong Zhang 
399b54502bSHong Zhang   PetscFunctionBegin;
40298cc208SBarry Smith   if (!ilu->inplace) {ierr = MatDestroy(&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);}
416bf464f9SBarry Smith   if (ilu->row && ilu->col && ilu->row != ilu->col) {ierr = ISDestroy(&ilu->row);CHKERRQ(ierr);}
42fcfd50ebSBarry Smith   ierr = ISDestroy(&ilu->col);CHKERRQ(ierr);
439b54502bSHong Zhang   PetscFunctionReturn(0);
449b54502bSHong Zhang }
459b54502bSHong Zhang 
469b54502bSHong Zhang #undef __FUNCT__
47b7c853c4SBarry Smith #define __FUNCT__ "PCFactorSetDropTolerance_ILU"
487087cfbeSBarry Smith PetscErrorCode  PCFactorSetDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
499b54502bSHong Zhang {
50075768bcSBarry Smith   PC_ILU *ilu = (PC_ILU*)pc->data;
519b54502bSHong Zhang 
529b54502bSHong Zhang   PetscFunctionBegin;
534c9036c7SBarry Smith   if (pc->setupcalled && (((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
54ce94432eSBarry Smith     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot change drop tolerance after using PC");
559b54502bSHong Zhang   }
56075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dt      = dt;
57075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcol   = dtcol;
58075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcount = dtcount;
594c9036c7SBarry Smith   ((PC_Factor*)ilu)->info.usedt   = 1.0;
609b54502bSHong Zhang   PetscFunctionReturn(0);
619b54502bSHong Zhang }
629b54502bSHong Zhang 
639b54502bSHong Zhang #undef __FUNCT__
642401956bSBarry Smith #define __FUNCT__ "PCFactorSetReuseOrdering_ILU"
657087cfbeSBarry Smith PetscErrorCode  PCFactorSetReuseOrdering_ILU(PC pc,PetscBool flag)
669b54502bSHong Zhang {
67075768bcSBarry Smith   PC_ILU *ilu = (PC_ILU*)pc->data;
689b54502bSHong Zhang 
699b54502bSHong Zhang   PetscFunctionBegin;
709b54502bSHong Zhang   ilu->reuseordering = flag;
719b54502bSHong Zhang   PetscFunctionReturn(0);
729b54502bSHong Zhang }
739b54502bSHong Zhang 
749b54502bSHong Zhang #undef __FUNCT__
752401956bSBarry Smith #define __FUNCT__ "PCFactorSetUseInPlace_ILU"
768e37d05fSBarry Smith PetscErrorCode  PCFactorSetUseInPlace_ILU(PC pc,PetscBool flg)
779b54502bSHong Zhang {
78075768bcSBarry Smith   PC_ILU *dir = (PC_ILU*)pc->data;
799b54502bSHong Zhang 
809b54502bSHong Zhang   PetscFunctionBegin;
818e37d05fSBarry Smith   dir->inplace = flg;
828e37d05fSBarry Smith   PetscFunctionReturn(0);
838e37d05fSBarry Smith }
848e37d05fSBarry Smith 
858e37d05fSBarry Smith #undef __FUNCT__
868e37d05fSBarry Smith #define __FUNCT__ "PCFactorGetUseInPlace_ILU"
878e37d05fSBarry Smith PetscErrorCode  PCFactorGetUseInPlace_ILU(PC pc,PetscBool *flg)
888e37d05fSBarry Smith {
898e37d05fSBarry Smith   PC_ILU *dir = (PC_ILU*)pc->data;
908e37d05fSBarry Smith 
918e37d05fSBarry Smith   PetscFunctionBegin;
928e37d05fSBarry Smith   *flg = dir->inplace;
939b54502bSHong Zhang   PetscFunctionReturn(0);
949b54502bSHong Zhang }
959b54502bSHong Zhang 
969b54502bSHong Zhang #undef __FUNCT__
979b54502bSHong Zhang #define __FUNCT__ "PCSetFromOptions_ILU"
984416b707SBarry Smith static PetscErrorCode PCSetFromOptions_ILU(PetscOptionItems *PetscOptionsObject,PC pc)
999b54502bSHong Zhang {
1009b54502bSHong Zhang   PetscErrorCode ierr;
10178fc6b22SHong Zhang   PetscInt       itmp;
1028afaa268SBarry Smith   PetscBool      flg,set;
1039b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
1049b54502bSHong Zhang   PetscReal      tol;
1059b54502bSHong Zhang 
1069b54502bSHong Zhang   PetscFunctionBegin;
107e55864a3SBarry Smith   ierr = PetscOptionsHead(PetscOptionsObject,"ILU Options");CHKERRQ(ierr);
108e55864a3SBarry Smith   ierr = PCSetFromOptions_Factor(PetscOptionsObject,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 
1138afaa268SBarry Smith   ierr = PetscOptionsBool("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",((PC_Factor*)ilu)->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr);
1148afaa268SBarry Smith   if (set) ((PC_Factor*)ilu)->info.diagonal_fill = (PetscReal) flg;
1152401956bSBarry Smith   ierr = PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);CHKERRQ(ierr);
1169b54502bSHong Zhang   if (flg) {
1179b54502bSHong Zhang     tol  = PETSC_DECIDE;
1182401956bSBarry Smith     ierr = PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);CHKERRQ(ierr);
1192401956bSBarry Smith     ierr = PCFactorReorderForNonzeroDiagonal(pc,tol);CHKERRQ(ierr);
1209b54502bSHong Zhang   }
1219b54502bSHong Zhang 
1229b54502bSHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
1239b54502bSHong Zhang   PetscFunctionReturn(0);
1249b54502bSHong Zhang }
1259b54502bSHong Zhang 
1269b54502bSHong Zhang #undef __FUNCT__
1279b54502bSHong Zhang #define __FUNCT__ "PCView_ILU"
1289b54502bSHong Zhang static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
1299b54502bSHong Zhang {
1309b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
1319b54502bSHong Zhang   PetscErrorCode ierr;
132ace3abfcSBarry Smith   PetscBool      iascii;
1339b54502bSHong Zhang 
1349b54502bSHong Zhang   PetscFunctionBegin;
135251f4c67SDmitry Karpeev   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1369b54502bSHong Zhang   if (iascii) {
137914a5d51SHong Zhang     if (ilu->inplace) {
138914a5d51SHong Zhang       ierr = PetscViewerASCIIPrintf(viewer,"  ILU: in-place factorization\n");CHKERRQ(ierr);
1399b54502bSHong Zhang     } else {
140914a5d51SHong Zhang       ierr = PetscViewerASCIIPrintf(viewer,"  ILU: out-of-place factorization\n");CHKERRQ(ierr);
141145b9266SHong Zhang     }
142145b9266SHong Zhang 
143914a5d51SHong Zhang     if (ilu->reusefill)     {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: Reusing fill from past factorization\n");CHKERRQ(ierr);}
144914a5d51SHong Zhang     if (ilu->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: Reusing reordering from past factorization\n");CHKERRQ(ierr);}
1459b54502bSHong Zhang   }
146914a5d51SHong Zhang   ierr = PCView_Factor(pc,viewer);CHKERRQ(ierr);
1479b54502bSHong Zhang   PetscFunctionReturn(0);
1489b54502bSHong Zhang }
1499b54502bSHong Zhang 
1509b54502bSHong Zhang #undef __FUNCT__
1519b54502bSHong Zhang #define __FUNCT__ "PCSetUp_ILU"
1529b54502bSHong Zhang static PetscErrorCode PCSetUp_ILU(PC pc)
1539b54502bSHong Zhang {
1549b54502bSHong Zhang   PetscErrorCode         ierr;
1559b54502bSHong Zhang   PC_ILU                 *ilu = (PC_ILU*)pc->data;
156f3a39becSBarry Smith   MatInfo                info;
157ace3abfcSBarry Smith   PetscBool              flg;
15800c67f3bSHong Zhang   const MatSolverPackage stype;
159*00e125f8SBarry Smith   MatFactorError         err;
1609b54502bSHong Zhang 
1619b54502bSHong Zhang   PetscFunctionBegin;
16292927226SBarry Smith   /* ugly hack to change default, since it is not support by some matrix types */
16392927226SBarry Smith   if (((PC_Factor*)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
16422d28d08SBarry Smith     ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);CHKERRQ(ierr);
16592927226SBarry Smith     if (!flg) {
16622d28d08SBarry Smith       ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&flg);CHKERRQ(ierr);
16792927226SBarry Smith       if (!flg) {
16892927226SBarry Smith         ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
169955c1f14SBarry Smith         PetscInfo(pc,"Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO\n");CHKERRQ(ierr);
17092927226SBarry Smith       }
17192927226SBarry Smith     }
17292927226SBarry Smith   }
17392927226SBarry Smith 
17484d44b13SHong Zhang   ierr = MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);CHKERRQ(ierr);
1759b54502bSHong Zhang   if (ilu->inplace) {
1769b54502bSHong Zhang     if (!pc->setupcalled) {
1779b54502bSHong Zhang 
1789b54502bSHong Zhang       /* In-place factorization only makes sense with the natural ordering,
1799b54502bSHong Zhang          so we only need to get the ordering once, even if nonzero structure changes */
180075768bcSBarry Smith       ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
1813bb1ff40SBarry Smith       if (ilu->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);CHKERRQ(ierr);}
1823bb1ff40SBarry Smith       if (ilu->col) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);CHKERRQ(ierr);}
1839b54502bSHong Zhang     }
1849b54502bSHong Zhang 
1859b54502bSHong Zhang     /* In place ILU only makes sense with fill factor of 1.0 because
1869b54502bSHong Zhang        cannot have levels of fill */
187075768bcSBarry Smith     ((PC_Factor*)ilu)->info.fill          = 1.0;
18875567043SBarry Smith     ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
1892fa5cd67SKarl Rupp 
190365a8a9eSBarry Smith     ierr = MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);CHKERRQ(ierr);
191*00e125f8SBarry Smith     ierr = MatFactorGetError(pc->pmat,&err);CHKERRQ(ierr);
192*00e125f8SBarry Smith     if (err) { /* Factor() fails */
193*00e125f8SBarry Smith       pc->failedreason = (PCFailedReason)err;
1946baea169SHong Zhang       PetscFunctionReturn(0);
1956baea169SHong Zhang     }
1966baea169SHong Zhang 
197075768bcSBarry Smith     ((PC_Factor*)ilu)->fact = pc->pmat;
198b33856dcSBarry Smith     /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
199b33856dcSBarry Smith     ierr = PetscObjectStateGet((PetscObject)pc->pmat,&pc->matstate);CHKERRQ(ierr);
2009b54502bSHong Zhang   } else {
2019b54502bSHong Zhang     if (!pc->setupcalled) {
2029b54502bSHong Zhang       /* first time in so compute reordering and symbolic factorization */
203075768bcSBarry Smith       ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
2043bb1ff40SBarry Smith       if (ilu->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);CHKERRQ(ierr);}
2053bb1ff40SBarry Smith       if (ilu->col) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);CHKERRQ(ierr);}
2069b54502bSHong Zhang       /*  Remove zeros along diagonal?     */
2079b54502bSHong Zhang       if (ilu->nonzerosalongdiagonal) {
2089b54502bSHong Zhang         ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr);
2099b54502bSHong Zhang       }
21086daa42cSBarry Smith       if (!((PC_Factor*)ilu)->fact) {
21111bfe483SHong Zhang         ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
212a1f19f5aSHong Zhang       }
213075768bcSBarry Smith       ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
214075768bcSBarry Smith       ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
2152fa5cd67SKarl Rupp 
216f3a39becSBarry Smith       ilu->actualfill = info.fill_ratio_needed;
2172fa5cd67SKarl Rupp 
2183bb1ff40SBarry Smith       ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
2199b54502bSHong Zhang     } else if (pc->flag != SAME_NONZERO_PATTERN) {
2209b54502bSHong Zhang       if (!ilu->reuseordering) {
2219b54502bSHong Zhang         /* compute a new ordering for the ILU */
222fcfd50ebSBarry Smith         ierr = ISDestroy(&ilu->row);CHKERRQ(ierr);
223fcfd50ebSBarry Smith         ierr = ISDestroy(&ilu->col);CHKERRQ(ierr);
224075768bcSBarry Smith         ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
2253bb1ff40SBarry Smith         if (ilu->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);CHKERRQ(ierr);}
2263bb1ff40SBarry Smith         if (ilu->col) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);CHKERRQ(ierr);}
2279b54502bSHong Zhang         /*  Remove zeros along diagonal?     */
2289b54502bSHong Zhang         if (ilu->nonzerosalongdiagonal) {
2299b54502bSHong Zhang           ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr);
2309b54502bSHong Zhang         }
2319b54502bSHong Zhang       }
2326bf464f9SBarry Smith       ierr = MatDestroy(&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
23314798fb4SJed Brown       ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
234075768bcSBarry Smith       ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
235075768bcSBarry Smith       ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
2362fa5cd67SKarl Rupp 
237f3a39becSBarry Smith       ilu->actualfill = info.fill_ratio_needed;
2382fa5cd67SKarl Rupp 
2393bb1ff40SBarry Smith       ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
2409b54502bSHong Zhang     }
241*00e125f8SBarry Smith     ierr = MatFactorGetError(((PC_Factor*)ilu)->fact,&err);CHKERRQ(ierr);
242*00e125f8SBarry Smith     if (err) { /* FactorSymbolic() fails */
243*00e125f8SBarry Smith       pc->failedreason = (PCFailedReason)err;
2446baea169SHong Zhang       PetscFunctionReturn(0);
2456baea169SHong Zhang     }
2466baea169SHong Zhang 
247075768bcSBarry Smith     ierr = MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
248*00e125f8SBarry Smith     ierr = MatFactorGetError(((PC_Factor*)ilu)->fact,&err);CHKERRQ(ierr);
249*00e125f8SBarry Smith     if (err) { /* FactorNumeric() fails */
250*00e125f8SBarry Smith       pc->failedreason = (PCFailedReason)err;
2516baea169SHong Zhang     }
2529b54502bSHong Zhang   }
25300c67f3bSHong Zhang 
25400c67f3bSHong Zhang   ierr = PCFactorGetMatSolverPackage(pc,&stype);CHKERRQ(ierr);
25500c67f3bSHong Zhang   if (!stype) {
256*00e125f8SBarry Smith     const MatSolverPackage solverpackage;
257*00e125f8SBarry Smith     ierr = MatFactorGetSolverPackage(((PC_Factor*)ilu)->fact,&solverpackage);CHKERRQ(ierr);
258*00e125f8SBarry Smith     ierr = PCFactorSetMatSolverPackage(pc,solverpackage);CHKERRQ(ierr);
25900c67f3bSHong Zhang   }
2609b54502bSHong Zhang   PetscFunctionReturn(0);
2619b54502bSHong Zhang }
2629b54502bSHong Zhang 
2639b54502bSHong Zhang #undef __FUNCT__
2649b54502bSHong Zhang #define __FUNCT__ "PCDestroy_ILU"
2659b54502bSHong Zhang static PetscErrorCode PCDestroy_ILU(PC pc)
2669b54502bSHong Zhang {
2679b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2689b54502bSHong Zhang   PetscErrorCode ierr;
2699b54502bSHong Zhang 
2709b54502bSHong Zhang   PetscFunctionBegin;
271574deadeSBarry Smith   ierr = PCReset_ILU(pc);CHKERRQ(ierr);
272503cfb0cSBarry Smith   ierr = PetscFree(((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr);
273503cfb0cSBarry Smith   ierr = PetscFree(((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
274c31cb41cSBarry Smith   ierr = PetscFree(pc->data);CHKERRQ(ierr);
2759b54502bSHong Zhang   PetscFunctionReturn(0);
2769b54502bSHong Zhang }
2779b54502bSHong Zhang 
2789b54502bSHong Zhang #undef __FUNCT__
2799b54502bSHong Zhang #define __FUNCT__ "PCApply_ILU"
2809b54502bSHong Zhang static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
2819b54502bSHong Zhang {
2829b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2839b54502bSHong Zhang   PetscErrorCode ierr;
2849b54502bSHong Zhang 
2859b54502bSHong Zhang   PetscFunctionBegin;
286075768bcSBarry Smith   ierr = MatSolve(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr);
2879b54502bSHong Zhang   PetscFunctionReturn(0);
2889b54502bSHong Zhang }
2899b54502bSHong Zhang 
2909b54502bSHong Zhang #undef __FUNCT__
2919b54502bSHong Zhang #define __FUNCT__ "PCApplyTranspose_ILU"
2929b54502bSHong Zhang static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
2939b54502bSHong Zhang {
2949b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2959b54502bSHong Zhang   PetscErrorCode ierr;
2969b54502bSHong Zhang 
2979b54502bSHong Zhang   PetscFunctionBegin;
298075768bcSBarry Smith   ierr = MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr);
2999b54502bSHong Zhang   PetscFunctionReturn(0);
3009b54502bSHong Zhang }
3019b54502bSHong Zhang 
302f0b9ad6cSBarry Smith #undef __FUNCT__
303f0b9ad6cSBarry Smith #define __FUNCT__ "PCApplySymmetricLeft_ILU"
304f0b9ad6cSBarry Smith static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
305f0b9ad6cSBarry Smith {
306f0b9ad6cSBarry Smith   PetscErrorCode ierr;
307f0b9ad6cSBarry Smith   PC_ILU         *icc = (PC_ILU*)pc->data;
308f0b9ad6cSBarry Smith 
309f0b9ad6cSBarry Smith   PetscFunctionBegin;
310f0b9ad6cSBarry Smith   ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr);
311f0b9ad6cSBarry Smith   PetscFunctionReturn(0);
312f0b9ad6cSBarry Smith }
313f0b9ad6cSBarry Smith 
314f0b9ad6cSBarry Smith #undef __FUNCT__
315f0b9ad6cSBarry Smith #define __FUNCT__ "PCApplySymmetricRight_ILU"
316f0b9ad6cSBarry Smith static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
317f0b9ad6cSBarry Smith {
318f0b9ad6cSBarry Smith   PetscErrorCode ierr;
319f0b9ad6cSBarry Smith   PC_ILU         *icc = (PC_ILU*)pc->data;
320f0b9ad6cSBarry Smith 
321f0b9ad6cSBarry Smith   PetscFunctionBegin;
322f0b9ad6cSBarry Smith   ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr);
323f0b9ad6cSBarry Smith   PetscFunctionReturn(0);
324f0b9ad6cSBarry Smith }
325f0b9ad6cSBarry Smith 
3269b54502bSHong Zhang /*MC
3279b54502bSHong Zhang      PCILU - Incomplete factorization preconditioners.
3289b54502bSHong Zhang 
3299b54502bSHong Zhang    Options Database Keys:
3302401956bSBarry Smith +  -pc_factor_levels <k> - number of levels of fill for ILU(k)
3312401956bSBarry Smith .  -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
3329b54502bSHong Zhang                       its factorization (overwrites original matrix)
3332401956bSBarry Smith .  -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
3342401956bSBarry Smith .  -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
33555ba2a51SBarry Smith .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
3362401956bSBarry Smith .  -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
3379b54502bSHong Zhang                                    this decreases the chance of getting a zero pivot
3382401956bSBarry Smith .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
339967c93d3SBarry Smith -  -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
3409b54502bSHong Zhang                              than 1 the diagonal blocks are factored with partial pivoting (this increases the
3419b54502bSHong Zhang                              stability of the ILU factorization
3429b54502bSHong Zhang 
3439b54502bSHong Zhang    Level: beginner
3449b54502bSHong Zhang 
3459b54502bSHong Zhang   Concepts: incomplete factorization
3469b54502bSHong Zhang 
347d9ba8547SBarry Smith    Notes: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)
3489b54502bSHong Zhang 
3499b54502bSHong Zhang           For BAIJ matrices this implements a point block ILU
3509b54502bSHong Zhang 
351f0b9ad6cSBarry Smith           The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
352f0b9ad6cSBarry Smith           even when the matrix is not symmetric since the U stores the diagonals of the factorization.
353f0b9ad6cSBarry Smith 
354cd0a26f6SPaul Mullowney           If you are using MATSEQAIJCUSPARSE matrices (or MATMPIAIJCUSPARESE matrices with block Jacobi), factorization
355cd0a26f6SPaul Mullowney           is never done on the GPU).
356ac793be5SBarry Smith 
357c582cd25SBarry Smith    References:
35896a0c994SBarry Smith +  1. - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
35996a0c994SBarry Smith    self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968.
36096a0c994SBarry Smith .  2. -  T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961.
36196a0c994SBarry Smith -  3. -  TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
36296a0c994SBarry Smith       Chapter in Parallel Numerical
363c582cd25SBarry Smith       Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
36496a0c994SBarry Smith       Science and Engineering, Kluwer.
365c582cd25SBarry Smith 
366c582cd25SBarry Smith 
3679b54502bSHong Zhang .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
368145b9266SHong Zhang            PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
369b7c853c4SBarry Smith            PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
37092e9c092SBarry Smith            PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(),
37192e9c092SBarry Smith            PCFactorGetAllowDiagonalFill(), PCFactorGetUseInPlace()
3729b54502bSHong Zhang 
3739b54502bSHong Zhang M*/
3749b54502bSHong Zhang 
3759b54502bSHong Zhang #undef __FUNCT__
3769b54502bSHong Zhang #define __FUNCT__ "PCCreate_ILU"
3778cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
3789b54502bSHong Zhang {
3799b54502bSHong Zhang   PetscErrorCode ierr;
3809b54502bSHong Zhang   PC_ILU         *ilu;
3819b54502bSHong Zhang 
3829b54502bSHong Zhang   PetscFunctionBegin;
383b00a9115SJed Brown   ierr = PetscNewLog(pc,&ilu);CHKERRQ(ierr);
3849b54502bSHong Zhang 
385075768bcSBarry Smith   ((PC_Factor*)ilu)->fact               = 0;
386075768bcSBarry Smith   ierr                                  = MatFactorInfoInitialize(&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
387879e8a4dSBarry Smith   ((PC_Factor*)ilu)->factortype         = MAT_FACTOR_ILU;
38875567043SBarry Smith   ((PC_Factor*)ilu)->info.levels        = 0.;
389075768bcSBarry Smith   ((PC_Factor*)ilu)->info.fill          = 1.0;
3909b54502bSHong Zhang   ilu->col                              = 0;
3919b54502bSHong Zhang   ilu->row                              = 0;
3929b54502bSHong Zhang   ilu->inplace                          = PETSC_FALSE;
39319fd82e9SBarry Smith   ierr                                  = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
3949b54502bSHong Zhang   ilu->reuseordering                    = PETSC_FALSE;
395075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dt            = PETSC_DEFAULT;
396075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcount       = PETSC_DEFAULT;
397075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcol         = PETSC_DEFAULT;
398c378d667SBarry Smith   ((PC_Factor*)ilu)->info.shifttype     = (PetscReal)MAT_SHIFT_NONE;
3990ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.shiftamount   = 100.0*PETSC_MACHINE_EPSILON;
4000ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.zeropivot     = 100.0*PETSC_MACHINE_EPSILON;
401075768bcSBarry Smith   ((PC_Factor*)ilu)->info.pivotinblocks = 1.0;
4029b54502bSHong Zhang   ilu->reusefill                        = PETSC_FALSE;
4030ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
4049b54502bSHong Zhang   pc->data                              = (void*)ilu;
4059b54502bSHong Zhang 
406574deadeSBarry Smith   pc->ops->reset               = PCReset_ILU;
4079b54502bSHong Zhang   pc->ops->destroy             = PCDestroy_ILU;
4089b54502bSHong Zhang   pc->ops->apply               = PCApply_ILU;
4099b54502bSHong Zhang   pc->ops->applytranspose      = PCApplyTranspose_ILU;
4109b54502bSHong Zhang   pc->ops->setup               = PCSetUp_ILU;
4119b54502bSHong Zhang   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
41285317021SBarry Smith   pc->ops->getfactoredmatrix   = PCFactorGetMatrix_Factor;
4139b54502bSHong Zhang   pc->ops->view                = PCView_ILU;
414f0b9ad6cSBarry Smith   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ILU;
415f0b9ad6cSBarry Smith   pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
4169b54502bSHong Zhang   pc->ops->applyrichardson     = 0;
4179b54502bSHong Zhang 
418bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);CHKERRQ(ierr);
419bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);CHKERRQ(ierr);
420bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);CHKERRQ(ierr);
421bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr);
422bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr);
423bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr);
424bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);CHKERRQ(ierr);
425bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);CHKERRQ(ierr);
426bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr);
427bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_ILU);CHKERRQ(ierr);
428bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_ILU);CHKERRQ(ierr);
429bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C",PCFactorSetLevels_Factor);CHKERRQ(ierr);
4302591b318SToby Isaac   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetLevels_C",PCFactorGetLevels_Factor);CHKERRQ(ierr);
431bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_ILU);CHKERRQ(ierr);
4328e37d05fSBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_ILU);CHKERRQ(ierr);
433bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C",PCFactorSetAllowDiagonalFill_Factor);CHKERRQ(ierr);
43492e9c092SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetAllowDiagonalFill_C",PCFactorGetAllowDiagonalFill_Factor);CHKERRQ(ierr);
435bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetPivotInBlocks_C",PCFactorSetPivotInBlocks_Factor);CHKERRQ(ierr);
436bdf89e91SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_ILU);CHKERRQ(ierr);
4379b54502bSHong Zhang   PetscFunctionReturn(0);
4389b54502bSHong Zhang }
439