xref: /petsc/src/ksp/pc/impls/factor/ilu/ilu.c (revision 8cc058d9cd56c1ccb3be12a47760ddfc446aaffc)
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"
767087cfbeSBarry Smith PetscErrorCode  PCFactorSetUseInPlace_ILU(PC pc)
779b54502bSHong Zhang {
78075768bcSBarry Smith   PC_ILU *dir = (PC_ILU*)pc->data;
799b54502bSHong Zhang 
809b54502bSHong Zhang   PetscFunctionBegin;
819b54502bSHong Zhang   dir->inplace = PETSC_TRUE;
829b54502bSHong Zhang   PetscFunctionReturn(0);
839b54502bSHong Zhang }
849b54502bSHong Zhang 
859b54502bSHong Zhang #undef __FUNCT__
869b54502bSHong Zhang #define __FUNCT__ "PCSetFromOptions_ILU"
879b54502bSHong Zhang static PetscErrorCode PCSetFromOptions_ILU(PC pc)
889b54502bSHong Zhang {
899b54502bSHong Zhang   PetscErrorCode ierr;
9078fc6b22SHong Zhang   PetscInt       itmp;
91ace3abfcSBarry Smith   PetscBool      flg;
929b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
939b54502bSHong Zhang   PetscReal      tol;
942fa5cd67SKarl Rupp   /* PetscReal      dt[3]; */
959b54502bSHong Zhang 
969b54502bSHong Zhang   PetscFunctionBegin;
979b54502bSHong Zhang   ierr = PetscOptionsHead("ILU Options");CHKERRQ(ierr);
988ff23777SHong Zhang   ierr = PCSetFromOptions_Factor(pc);CHKERRQ(ierr);
998ff23777SHong Zhang 
100075768bcSBarry Smith   ierr = PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);CHKERRQ(ierr);
101075768bcSBarry Smith   if (flg) ((PC_Factor*)ilu)->info.levels = itmp;
1022fa5cd67SKarl Rupp 
10390d69ab7SBarry Smith   flg  = PETSC_FALSE;
1040298fd71SBarry Smith   ierr = PetscOptionsBool("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",flg,&flg,NULL);CHKERRQ(ierr);
105075768bcSBarry Smith   ((PC_Factor*)ilu)->info.diagonal_fill = (double) flg;
1065a586d82SBarry Smith   /*
107075768bcSBarry Smith   dt[0] = ((PC_Factor*)ilu)->info.dt;
108075768bcSBarry Smith   dt[1] = ((PC_Factor*)ilu)->info.dtcol;
109075768bcSBarry Smith   dt[2] = ((PC_Factor*)ilu)->info.dtcount;
1105a586d82SBarry Smith 
11178fc6b22SHong Zhang   PetscInt       dtmax = 3;
11278fc6b22SHong Zhang   ierr = PetscOptionsRealArray("-pc_factor_drop_tolerance,","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);CHKERRQ(ierr);
1139b54502bSHong Zhang   if (flg) {
114b7c853c4SBarry Smith     ierr = PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);CHKERRQ(ierr);
1159b54502bSHong Zhang   }
11678fc6b22SHong Zhang   */
1172401956bSBarry Smith   ierr = PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);CHKERRQ(ierr);
1189b54502bSHong Zhang   if (flg) {
1199b54502bSHong Zhang     tol  = PETSC_DECIDE;
1202401956bSBarry Smith     ierr = PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);CHKERRQ(ierr);
1212401956bSBarry Smith     ierr = PCFactorReorderForNonzeroDiagonal(pc,tol);CHKERRQ(ierr);
1229b54502bSHong Zhang   }
1239b54502bSHong Zhang 
1249b54502bSHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
1259b54502bSHong Zhang   PetscFunctionReturn(0);
1269b54502bSHong Zhang }
1279b54502bSHong Zhang 
1289b54502bSHong Zhang #undef __FUNCT__
1299b54502bSHong Zhang #define __FUNCT__ "PCView_ILU"
1309b54502bSHong Zhang static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
1319b54502bSHong Zhang {
1329b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
1339b54502bSHong Zhang   PetscErrorCode ierr;
134ace3abfcSBarry Smith   PetscBool      iascii;
1359b54502bSHong Zhang 
1369b54502bSHong Zhang   PetscFunctionBegin;
137251f4c67SDmitry Karpeev   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1389b54502bSHong Zhang   if (iascii) {
139914a5d51SHong Zhang     if (ilu->inplace) {
140914a5d51SHong Zhang       ierr = PetscViewerASCIIPrintf(viewer,"  ILU: in-place factorization\n");CHKERRQ(ierr);
1419b54502bSHong Zhang     } else {
142914a5d51SHong Zhang       ierr = PetscViewerASCIIPrintf(viewer,"  ILU: out-of-place factorization\n");CHKERRQ(ierr);
143145b9266SHong Zhang     }
144145b9266SHong Zhang 
145914a5d51SHong Zhang     if (ilu->reusefill)     {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: Reusing fill from past factorization\n");CHKERRQ(ierr);}
146914a5d51SHong Zhang     if (ilu->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer,"  ILU: Reusing reordering from past factorization\n");CHKERRQ(ierr);}
1479b54502bSHong Zhang   }
148914a5d51SHong Zhang   ierr = PCView_Factor(pc,viewer);CHKERRQ(ierr);
1499b54502bSHong Zhang   PetscFunctionReturn(0);
1509b54502bSHong Zhang }
1519b54502bSHong Zhang 
1529b54502bSHong Zhang #undef __FUNCT__
1539b54502bSHong Zhang #define __FUNCT__ "PCSetUp_ILU"
1549b54502bSHong Zhang static PetscErrorCode PCSetUp_ILU(PC pc)
1559b54502bSHong Zhang {
1569b54502bSHong Zhang   PetscErrorCode ierr;
1579b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
158f3a39becSBarry Smith   MatInfo        info;
159ace3abfcSBarry Smith   PetscBool      flg;
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;
16922d28d08SBarry Smith         PetscInfo(pc,"Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO");CHKERRQ(ierr);
17092927226SBarry Smith       }
17192927226SBarry Smith     }
17292927226SBarry Smith   }
17392927226SBarry Smith 
1749b54502bSHong Zhang   if (ilu->inplace) {
175075768bcSBarry Smith     CHKMEMQ;
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);
181efee365bSSatish Balay       if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
182efee365bSSatish Balay       if (ilu->col) {ierr = PetscLogObjectParent(pc,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 
190fe97e370SBarry Smith     ierr = MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);CHKMEMQ;
191075768bcSBarry Smith     ((PC_Factor*)ilu)->fact = pc->pmat;
1929b54502bSHong Zhang   } else {
1939b54502bSHong Zhang     if (!pc->setupcalled) {
1949b54502bSHong Zhang       /* first time in so compute reordering and symbolic factorization */
195075768bcSBarry Smith       ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
196efee365bSSatish Balay       if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
197efee365bSSatish Balay       if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
1989b54502bSHong Zhang       /*  Remove zeros along diagonal?     */
1999b54502bSHong Zhang       if (ilu->nonzerosalongdiagonal) {
2009b54502bSHong Zhang         ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr);
2019b54502bSHong Zhang       }
20286daa42cSBarry Smith       if (!((PC_Factor*)ilu)->fact) {
20311bfe483SHong Zhang         ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
204a1f19f5aSHong Zhang       }
205075768bcSBarry Smith       ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
206075768bcSBarry Smith       ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
2072fa5cd67SKarl Rupp 
208f3a39becSBarry Smith       ilu->actualfill = info.fill_ratio_needed;
2092fa5cd67SKarl Rupp 
210075768bcSBarry Smith       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
2119b54502bSHong Zhang     } else if (pc->flag != SAME_NONZERO_PATTERN) {
2129b54502bSHong Zhang       if (!ilu->reuseordering) {
2139b54502bSHong Zhang         /* compute a new ordering for the ILU */
214fcfd50ebSBarry Smith         ierr = ISDestroy(&ilu->row);CHKERRQ(ierr);
215fcfd50ebSBarry Smith         ierr = ISDestroy(&ilu->col);CHKERRQ(ierr);
216075768bcSBarry Smith         ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
217efee365bSSatish Balay         if (ilu->row) {ierr = PetscLogObjectParent(pc,ilu->row);CHKERRQ(ierr);}
218efee365bSSatish Balay         if (ilu->col) {ierr = PetscLogObjectParent(pc,ilu->col);CHKERRQ(ierr);}
2199b54502bSHong Zhang         /*  Remove zeros along diagonal?     */
2209b54502bSHong Zhang         if (ilu->nonzerosalongdiagonal) {
2219b54502bSHong Zhang           ierr = MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);CHKERRQ(ierr);
2229b54502bSHong Zhang         }
2239b54502bSHong Zhang       }
2246bf464f9SBarry Smith       ierr = MatDestroy(&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
22514798fb4SJed Brown       ierr = MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
226075768bcSBarry Smith       ierr = MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
227075768bcSBarry Smith       ierr = MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
2282fa5cd67SKarl Rupp 
229f3a39becSBarry Smith       ilu->actualfill = info.fill_ratio_needed;
2302fa5cd67SKarl Rupp 
231075768bcSBarry Smith       ierr = PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);CHKERRQ(ierr);
2329b54502bSHong Zhang     }
233075768bcSBarry Smith     CHKMEMQ;
234075768bcSBarry Smith     ierr = MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
235075768bcSBarry Smith     CHKMEMQ;
2369b54502bSHong Zhang   }
2379b54502bSHong Zhang   PetscFunctionReturn(0);
2389b54502bSHong Zhang }
2399b54502bSHong Zhang 
2409b54502bSHong Zhang #undef __FUNCT__
2419b54502bSHong Zhang #define __FUNCT__ "PCDestroy_ILU"
2429b54502bSHong Zhang static PetscErrorCode PCDestroy_ILU(PC pc)
2439b54502bSHong Zhang {
2449b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2459b54502bSHong Zhang   PetscErrorCode ierr;
2469b54502bSHong Zhang 
2479b54502bSHong Zhang   PetscFunctionBegin;
248574deadeSBarry Smith   ierr = PCReset_ILU(pc);CHKERRQ(ierr);
249503cfb0cSBarry Smith   ierr = PetscFree(((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr);
250503cfb0cSBarry Smith   ierr = PetscFree(((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
251c31cb41cSBarry Smith   ierr = PetscFree(pc->data);CHKERRQ(ierr);
2529b54502bSHong Zhang   PetscFunctionReturn(0);
2539b54502bSHong Zhang }
2549b54502bSHong Zhang 
2559b54502bSHong Zhang #undef __FUNCT__
2569b54502bSHong Zhang #define __FUNCT__ "PCApply_ILU"
2579b54502bSHong Zhang static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
2589b54502bSHong Zhang {
2599b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2609b54502bSHong Zhang   PetscErrorCode ierr;
2619b54502bSHong Zhang 
2629b54502bSHong Zhang   PetscFunctionBegin;
263075768bcSBarry Smith   ierr = MatSolve(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr);
2649b54502bSHong Zhang   PetscFunctionReturn(0);
2659b54502bSHong Zhang }
2669b54502bSHong Zhang 
2679b54502bSHong Zhang #undef __FUNCT__
2689b54502bSHong Zhang #define __FUNCT__ "PCApplyTranspose_ILU"
2699b54502bSHong Zhang static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
2709b54502bSHong Zhang {
2719b54502bSHong Zhang   PC_ILU         *ilu = (PC_ILU*)pc->data;
2729b54502bSHong Zhang   PetscErrorCode ierr;
2739b54502bSHong Zhang 
2749b54502bSHong Zhang   PetscFunctionBegin;
275075768bcSBarry Smith   ierr = MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);CHKERRQ(ierr);
2769b54502bSHong Zhang   PetscFunctionReturn(0);
2779b54502bSHong Zhang }
2789b54502bSHong Zhang 
279f0b9ad6cSBarry Smith #undef __FUNCT__
280f0b9ad6cSBarry Smith #define __FUNCT__ "PCApplySymmetricLeft_ILU"
281f0b9ad6cSBarry Smith static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
282f0b9ad6cSBarry Smith {
283f0b9ad6cSBarry Smith   PetscErrorCode ierr;
284f0b9ad6cSBarry Smith   PC_ILU         *icc = (PC_ILU*)pc->data;
285f0b9ad6cSBarry Smith 
286f0b9ad6cSBarry Smith   PetscFunctionBegin;
287f0b9ad6cSBarry Smith   ierr = MatForwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr);
288f0b9ad6cSBarry Smith   PetscFunctionReturn(0);
289f0b9ad6cSBarry Smith }
290f0b9ad6cSBarry Smith 
291f0b9ad6cSBarry Smith #undef __FUNCT__
292f0b9ad6cSBarry Smith #define __FUNCT__ "PCApplySymmetricRight_ILU"
293f0b9ad6cSBarry Smith static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
294f0b9ad6cSBarry Smith {
295f0b9ad6cSBarry Smith   PetscErrorCode ierr;
296f0b9ad6cSBarry Smith   PC_ILU         *icc = (PC_ILU*)pc->data;
297f0b9ad6cSBarry Smith 
298f0b9ad6cSBarry Smith   PetscFunctionBegin;
299f0b9ad6cSBarry Smith   ierr = MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);CHKERRQ(ierr);
300f0b9ad6cSBarry Smith   PetscFunctionReturn(0);
301f0b9ad6cSBarry Smith }
302f0b9ad6cSBarry Smith 
3039b54502bSHong Zhang /*MC
3049b54502bSHong Zhang      PCILU - Incomplete factorization preconditioners.
3059b54502bSHong Zhang 
3069b54502bSHong Zhang    Options Database Keys:
3072401956bSBarry Smith +  -pc_factor_levels <k> - number of levels of fill for ILU(k)
3082401956bSBarry Smith .  -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
3099b54502bSHong Zhang                       its factorization (overwrites original matrix)
3102401956bSBarry Smith .  -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
3112401956bSBarry Smith .  -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
31255ba2a51SBarry Smith .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
3132401956bSBarry Smith .  -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
3149b54502bSHong Zhang                                    this decreases the chance of getting a zero pivot
3152401956bSBarry Smith .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
3162401956bSBarry Smith .  -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
3179b54502bSHong Zhang                              than 1 the diagonal blocks are factored with partial pivoting (this increases the
3189b54502bSHong Zhang                              stability of the ILU factorization
3199b54502bSHong Zhang 
3209b54502bSHong Zhang    Level: beginner
3219b54502bSHong Zhang 
3229b54502bSHong Zhang   Concepts: incomplete factorization
3239b54502bSHong Zhang 
324d9ba8547SBarry Smith    Notes: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)
3259b54502bSHong Zhang 
3269b54502bSHong Zhang           For BAIJ matrices this implements a point block ILU
3279b54502bSHong Zhang 
328f0b9ad6cSBarry Smith           The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
329f0b9ad6cSBarry Smith           even when the matrix is not symmetric since the U stores the diagonals of the factorization.
330f0b9ad6cSBarry Smith 
331ac793be5SBarry Smith           If you are using MATSEQAIJCUSP matrices (or MATMPIAIJCUSP matrices with block Jacobi) you must have ./configured
332ac793be5SBarry Smith           PETSc with also --download-txpetscgpu to have the triangular solves performed on the GPU (factorization is never
333ac793be5SBarry Smith           done on the GPU).
334ac793be5SBarry Smith 
335c582cd25SBarry Smith    References:
336c582cd25SBarry Smith    T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
337c582cd25SBarry Smith    self-adjoint elliptic difference equations. SIAM J. Numer. Anal., 5:559--573, 1968.
338c582cd25SBarry Smith 
339c582cd25SBarry Smith    T.A. Oliphant. An implicit numerical method for solving two-dimensional time-dependent dif-
340c582cd25SBarry Smith    fusion problems. Quart. Appl. Math., 19:221--229, 1961.
341c582cd25SBarry Smith 
342c582cd25SBarry Smith    Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST
343c582cd25SBarry Smith       http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical
344c582cd25SBarry Smith       Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
345c582cd25SBarry Smith       Science and Engineering, Kluwer, pp. 167--202.
346c582cd25SBarry Smith 
347c582cd25SBarry Smith 
3489b54502bSHong Zhang .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
349145b9266SHong Zhang            PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
350b7c853c4SBarry Smith            PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
3518ff23777SHong Zhang            PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks()
3529b54502bSHong Zhang 
3539b54502bSHong Zhang M*/
3549b54502bSHong Zhang 
3559b54502bSHong Zhang #undef __FUNCT__
3569b54502bSHong Zhang #define __FUNCT__ "PCCreate_ILU"
357*8cc058d9SJed Brown PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
3589b54502bSHong Zhang {
3599b54502bSHong Zhang   PetscErrorCode ierr;
3609b54502bSHong Zhang   PC_ILU         *ilu;
3619b54502bSHong Zhang 
3629b54502bSHong Zhang   PetscFunctionBegin;
36338f2d2fdSLisandro Dalcin   ierr = PetscNewLog(pc,PC_ILU,&ilu);CHKERRQ(ierr);
3649b54502bSHong Zhang 
365075768bcSBarry Smith   ((PC_Factor*)ilu)->fact               = 0;
366075768bcSBarry Smith   ierr                                  = MatFactorInfoInitialize(&((PC_Factor*)ilu)->info);CHKERRQ(ierr);
367879e8a4dSBarry Smith   ((PC_Factor*)ilu)->factortype         = MAT_FACTOR_ILU;
36875567043SBarry Smith   ((PC_Factor*)ilu)->info.levels        = 0.;
369075768bcSBarry Smith   ((PC_Factor*)ilu)->info.fill          = 1.0;
3709b54502bSHong Zhang   ilu->col                              = 0;
3719b54502bSHong Zhang   ilu->row                              = 0;
3729b54502bSHong Zhang   ilu->inplace                          = PETSC_FALSE;
3732692d6eeSBarry Smith   ierr                                  = PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)ilu)->solvertype);CHKERRQ(ierr);
37419fd82e9SBarry Smith   ierr                                  = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)ilu)->ordering);CHKERRQ(ierr);
3759b54502bSHong Zhang   ilu->reuseordering                    = PETSC_FALSE;
376075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dt            = PETSC_DEFAULT;
377075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcount       = PETSC_DEFAULT;
378075768bcSBarry Smith   ((PC_Factor*)ilu)->info.dtcol         = PETSC_DEFAULT;
3799dbfc187SHong Zhang   ((PC_Factor*)ilu)->info.shifttype     = (PetscReal)MAT_SHIFT_INBLOCKS;
3800ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.shiftamount   = 100.0*PETSC_MACHINE_EPSILON;
3810ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.zeropivot     = 100.0*PETSC_MACHINE_EPSILON;
382075768bcSBarry Smith   ((PC_Factor*)ilu)->info.pivotinblocks = 1.0;
3839b54502bSHong Zhang   ilu->reusefill                        = PETSC_FALSE;
3840ed735ceSHong Zhang   ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
3859b54502bSHong Zhang   pc->data                              = (void*)ilu;
3869b54502bSHong Zhang 
387574deadeSBarry Smith   pc->ops->reset               = PCReset_ILU;
3889b54502bSHong Zhang   pc->ops->destroy             = PCDestroy_ILU;
3899b54502bSHong Zhang   pc->ops->apply               = PCApply_ILU;
3909b54502bSHong Zhang   pc->ops->applytranspose      = PCApplyTranspose_ILU;
3919b54502bSHong Zhang   pc->ops->setup               = PCSetUp_ILU;
3929b54502bSHong Zhang   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
39385317021SBarry Smith   pc->ops->getfactoredmatrix   = PCFactorGetMatrix_Factor;
3949b54502bSHong Zhang   pc->ops->view                = PCView_ILU;
395f0b9ad6cSBarry Smith   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ILU;
396f0b9ad6cSBarry Smith   pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
3979b54502bSHong Zhang   pc->ops->applyrichardson     = 0;
3989b54502bSHong Zhang 
39900de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",PCFactorSetZeroPivot_Factor);CHKERRQ(ierr);
40000de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor",PCFactorSetShiftType_Factor);CHKERRQ(ierr);
40100de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor",PCFactorSetShiftAmount_Factor);CHKERRQ(ierr);
40200de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr);
40300de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr);
40400de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor",PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr);
40500de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C","PCFactorSetDropTolerance_ILU",PCFactorSetDropTolerance_ILU);CHKERRQ(ierr);
40600de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",PCFactorSetFill_Factor);CHKERRQ(ierr);
40700de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr);
40800de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_ILU",PCFactorSetReuseOrdering_ILU);CHKERRQ(ierr);
40900de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_ILU",PCFactorSetReuseFill_ILU);CHKERRQ(ierr);
41000de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor",PCFactorSetLevels_Factor);CHKERRQ(ierr);
41100de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_ILU",PCFactorSetUseInPlace_ILU);CHKERRQ(ierr);
41200de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C","PCFactorSetAllowDiagonalFill_Factor",PCFactorSetAllowDiagonalFill_Factor);CHKERRQ(ierr);
41300de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetPivotInBlocks_C","PCFactorSetPivotInBlocks_Factor",PCFactorSetPivotInBlocks_Factor);CHKERRQ(ierr);
41400de8ff0SBarry Smith   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C","PCFactorReorderForNonzeroDiagonal_ILU",PCFactorReorderForNonzeroDiagonal_ILU);CHKERRQ(ierr);
4159b54502bSHong Zhang   PetscFunctionReturn(0);
4169b54502bSHong Zhang }
417