Actual source code: ilu.c

      1: #define PETSCKSP_DLL

      3: /*
      4:    Defines a ILU factorization preconditioner for any Mat implementation
      5: */
     6:  #include src/ksp/pc/pcimpl.h
     7:  #include src/ksp/pc/impls/factor/ilu/ilu.h
     8:  #include src/mat/matimpl.h

     10: /* ------------------------------------------------------------------------------------------*/
     14: PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetZeroPivot_ILU(PC pc,PetscReal z)
     15: {
     16:   PC_ILU *ilu;

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

     28: PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftNonzero_ILU(PC pc,PetscReal shift)
     29: {
     30:   PC_ILU *dir;

     33:   dir = (PC_ILU*)pc->data;
     34:   if (shift == (PetscReal) PETSC_DECIDE) {
     35:     dir->info.shiftnz = 1.e-12;
     36:   } else {
     37:     dir->info.shiftnz = shift;
     38:   }
     39:   return(0);
     40: }

     46: PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftPd_ILU(PC pc,PetscTruth shift)
     47: {
     48:   PC_ILU *dir;
     49:
     51:   dir = (PC_ILU*)pc->data;
     52:   dir->info.shiftpd = shift;
     53:   if (shift) dir->info.shift_fraction = 0.0;
     54:   return(0);
     55: }

     61: PetscErrorCode PETSCKSP_DLLEXPORT PCILUReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)
     62: {
     63:   PC_ILU *ilu = (PC_ILU*)pc->data;

     66:   ilu->nonzerosalongdiagonal = PETSC_TRUE;
     67:   if (z == PETSC_DECIDE) {
     68:     ilu->nonzerosalongdiagonaltol = 1.e-10;
     69:   } else {
     70:     ilu->nonzerosalongdiagonaltol = z;
     71:   }
     72:   return(0);
     73: }

     78: PetscErrorCode PCDestroy_ILU_Internal(PC pc)
     79: {
     80:   PC_ILU         *ilu = (PC_ILU*)pc->data;

     84:   if (!ilu->inplace && ilu->fact) {MatDestroy(ilu->fact);}
     85:   if (ilu->row && ilu->col && ilu->row != ilu->col) {ISDestroy(ilu->row);}
     86:   if (ilu->col) {ISDestroy(ilu->col);}
     87:   return(0);
     88: }

     93: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetUseDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
     94: {
     95:   PC_ILU         *ilu;

     99:   ilu = (PC_ILU*)pc->data;
    100:   if (pc->setupcalled && (!ilu->usedt || ilu->info.dt != dt || ilu->info.dtcol != dtcol || ilu->info.dtcount != dtcount)) {
    101:     pc->setupcalled   = 0;
    102:     PCDestroy_ILU_Internal(pc);
    103:   }
    104:   ilu->usedt        = PETSC_TRUE;
    105:   ilu->info.dt      = dt;
    106:   ilu->info.dtcol   = dtcol;
    107:   ilu->info.dtcount = dtcount;
    108:   ilu->info.fill    = PETSC_DEFAULT;
    109:   return(0);
    110: }

    116: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetFill_ILU(PC pc,PetscReal fill)
    117: {
    118:   PC_ILU *dir;

    121:   dir            = (PC_ILU*)pc->data;
    122:   dir->info.fill = fill;
    123:   return(0);
    124: }

    130: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetMatOrdering_ILU(PC pc,MatOrderingType ordering)
    131: {
    132:   PC_ILU         *dir = (PC_ILU*)pc->data;
    134:   PetscTruth     flg;
    135:
    137:   if (!pc->setupcalled) {
    138:      PetscStrfree(dir->ordering);
    139:      PetscStrallocpy(ordering,&dir->ordering);
    140:   } else {
    141:     PetscStrcmp(dir->ordering,ordering,&flg);
    142:     if (!flg) {
    143:       pc->setupcalled = 0;
    144:       PetscStrfree(dir->ordering);
    145:       PetscStrallocpy(ordering,&dir->ordering);
    146:       /* free the data structures, then create them again */
    147:       PCDestroy_ILU_Internal(pc);
    148:     }
    149:   }
    150:   return(0);
    151: }

    157: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetReuseOrdering_ILU(PC pc,PetscTruth flag)
    158: {
    159:   PC_ILU *ilu;

    162:   ilu                = (PC_ILU*)pc->data;
    163:   ilu->reuseordering = flag;
    164:   return(0);
    165: }

    171: PetscErrorCode PETSCKSP_DLLEXPORT PCILUDTSetReuseFill_ILUDT(PC pc,PetscTruth flag)
    172: {
    173:   PC_ILU *ilu;

    176:   ilu = (PC_ILU*)pc->data;
    177:   ilu->reusefill = flag;
    178:   if (flag) SETERRQ(PETSC_ERR_SUP,"Not yet supported");
    179:   return(0);
    180: }

    186: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetLevels_ILU(PC pc,PetscInt levels)
    187: {
    188:   PC_ILU         *ilu;

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

    194:   if (!pc->setupcalled) {
    195:     ilu->info.levels = levels;
    196:   } else if (ilu->usedt || ilu->info.levels != levels) {
    197:     ilu->info.levels = levels;
    198:     pc->setupcalled  = 0;
    199:     ilu->usedt       = PETSC_FALSE;
    200:     PCDestroy_ILU_Internal(pc);
    201:   }
    202:   return(0);
    203: }

    209: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetUseInPlace_ILU(PC pc)
    210: {
    211:   PC_ILU *dir;

    214:   dir          = (PC_ILU*)pc->data;
    215:   dir->inplace = PETSC_TRUE;
    216:   return(0);
    217: }

    223: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetAllowDiagonalFill_ILU(PC pc)
    224: {
    225:   PC_ILU *dir;

    228:   dir                 = (PC_ILU*)pc->data;
    229:   dir->info.diagonal_fill = 1;
    230:   return(0);
    231: }

    236: /*@
    237:    PCILUSetUseDropTolerance - The preconditioner will use an ILU
    238:    based on a drop tolerance.

    240:    Collective on PC

    242:    Input Parameters:
    243: +  pc - the preconditioner context
    244: .  dt - the drop tolerance, try from 1.e-10 to .1
    245: .  dtcol - tolerance for column pivot, good values [0.1 to 0.01]
    246: -  maxrowcount - the max number of nonzeros allowed in a row, best value
    247:                  depends on the number of nonzeros in row of original matrix

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

    252:    Level: intermediate

    254:     Notes:
    255:       This uses the iludt() code of Saad's SPARSKIT package

    257: .keywords: PC, levels, reordering, factorization, incomplete, ILU
    258: @*/
    259: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetUseDropTolerance(PC pc,PetscReal dt,PetscReal dtcol,PetscInt maxrowcount)
    260: {
    261:   PetscErrorCode ierr,(*f)(PC,PetscReal,PetscReal,PetscInt);

    265:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetUseDropTolerance_C",(void (**)(void))&f);
    266:   if (f) {
    267:     (*f)(pc,dt,dtcol,maxrowcount);
    268:   }
    269:   return(0);
    270: }

    274: /*@
    275:    PCILUSetFill - Indicate the amount of fill you expect in the factored matrix,
    276:    where fill = number nonzeros in factor/number nonzeros in original matrix.

    278:    Collective on PC

    280:    Input Parameters:
    281: +  pc - the preconditioner context
    282: -  fill - amount of expected fill

    284:    Options Database Key:
    285: $  -pc_ilu_fill <fill>

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

    293:    Level: intermediate

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

    297: .seealso: PCLUSetFill()
    298: @*/
    299: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetFill(PC pc,PetscReal fill)
    300: {
    301:   PetscErrorCode ierr,(*f)(PC,PetscReal);

    305:   if (fill < 1.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Fill factor cannot be less than 1.0");
    306:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetFill_C",(void (**)(void))&f);
    307:   if (f) {
    308:     (*f)(pc,fill);
    309:   }
    310:   return(0);
    311: }

    315: /*@C
    316:     PCILUSetMatOrdering - Sets the ordering routine (to reduce fill) to
    317:     be used in the ILU factorization.

    319:     Collective on PC

    321:     Input Parameters:
    322: +   pc - the preconditioner context
    323: -   ordering - the matrix ordering name, for example, MATORDERING_ND or MATORDERING_RCM

    325:     Options Database Key:
    326: .   -pc_ilu_mat_ordering_type <nd,rcm,...> - Sets ordering routine

    328:     Level: intermediate

    330:     Notes: natural ordering is used by default

    332: .seealso: PCLUSetMatOrdering()

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

    336: @*/
    337: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetMatOrdering(PC pc,MatOrderingType ordering)
    338: {
    339:   PetscErrorCode ierr,(*f)(PC,MatOrderingType);

    343:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetMatOrdering_C",(void (**)(void))&f);
    344:   if (f) {
    345:     (*f)(pc,ordering);
    346:   }
    347:   return(0);
    348: }

    352: /*@
    353:    PCILUSetReuseOrdering - When similar matrices are factored, this
    354:    causes the ordering computed in the first factor to be used for all
    355:    following factors; applies to both fill and drop tolerance ILUs.

    357:    Collective on PC

    359:    Input Parameters:
    360: +  pc - the preconditioner context
    361: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

    363:    Options Database Key:
    364: .  -pc_ilu_reuse_ordering - Activate PCILUSetReuseOrdering()

    366:    Level: intermediate

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

    370: .seealso: PCILUDTSetReuseFill(), PCLUSetReuseOrdering(), PCLUSetReuseFill()
    371: @*/
    372: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetReuseOrdering(PC pc,PetscTruth flag)
    373: {
    374:   PetscErrorCode ierr,(*f)(PC,PetscTruth);

    378:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetReuseOrdering_C",(void (**)(void))&f);
    379:   if (f) {
    380:     (*f)(pc,flag);
    381:   }
    382:   return(0);
    383: }

    387: /*@
    388:    PCILUDTSetReuseFill - When matrices with same nonzero structure are ILUDT factored,
    389:    this causes later ones to use the fill computed in the initial factorization.

    391:    Collective on PC

    393:    Input Parameters:
    394: +  pc - the preconditioner context
    395: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

    397:    Options Database Key:
    398: .  -pc_iludt_reuse_fill - Activates PCILUDTSetReuseFill()

    400:    Level: intermediate

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

    404: .seealso: PCILUSetReuseOrdering(), PCLUSetReuseOrdering(), PCLUSetReuseFill()
    405: @*/
    406: PetscErrorCode PETSCKSP_DLLEXPORT PCILUDTSetReuseFill(PC pc,PetscTruth flag)
    407: {
    408:   PetscErrorCode ierr,(*f)(PC,PetscTruth);

    412:   PetscObjectQueryFunction((PetscObject)pc,"PCILUDTSetReuseFill_C",(void (**)(void))&f);
    413:   if (f) {
    414:     (*f)(pc,flag);
    415:   }
    416:   return(0);
    417: }

    421: /*@
    422:    PCILUSetLevels - Sets the number of levels of fill to use.

    424:    Collective on PC

    426:    Input Parameters:
    427: +  pc - the preconditioner context
    428: -  levels - number of levels of fill

    430:    Options Database Key:
    431: .  -pc_ilu_levels <levels> - Sets fill level

    433:    Level: intermediate

    435: .keywords: PC, levels, fill, factorization, incomplete, ILU
    436: @*/
    437: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetLevels(PC pc,PetscInt levels)
    438: {
    439:   PetscErrorCode ierr,(*f)(PC,PetscInt);

    443:   if (levels < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"negative levels");
    444:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetLevels_C",(void (**)(void))&f);
    445:   if (f) {
    446:     (*f)(pc,levels);
    447:   }
    448:   return(0);
    449: }

    453: /*@
    454:    PCILUSetAllowDiagonalFill - Causes all diagonal matrix entries to be
    455:    treated as level 0 fill even if there is no non-zero location.

    457:    Collective on PC

    459:    Input Parameters:
    460: +  pc - the preconditioner context

    462:    Options Database Key:
    463: .  -pc_ilu_diagonal_fill

    465:    Notes:
    466:    Does not apply with 0 fill.

    468:    Level: intermediate

    470: .keywords: PC, levels, fill, factorization, incomplete, ILU
    471: @*/
    472: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetAllowDiagonalFill(PC pc)
    473: {
    474:   PetscErrorCode ierr,(*f)(PC);

    478:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetAllowDiagonalFill_C",(void (**)(void))&f);
    479:   if (f) {
    480:     (*f)(pc);
    481:   }
    482:   return(0);
    483: }

    487: /*@
    488:    PCILUSetUseInPlace - Tells the system to do an in-place incomplete factorization.
    489:    Collective on PC

    491:    Input Parameters:
    492: .  pc - the preconditioner context

    494:    Options Database Key:
    495: .  -pc_ilu_in_place - Activates in-place factorization

    497:    Notes:
    498:    PCILUSetUseInPlace() is intended for use with matrix-free variants of
    499:    Krylov methods, or when a different matrices are employed for the linear
    500:    system and preconditioner, or with ASM preconditioning.  Do NOT use
    501:    this option if the linear system
    502:    matrix also serves as the preconditioning matrix, since the factored
    503:    matrix would then overwrite the original matrix.

    505:    Only works well with ILU(0).

    507:    Level: intermediate

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

    511: .seealso:  PCLUSetUseInPlace()
    512: @*/
    513: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetUseInPlace(PC pc)
    514: {
    515:   PetscErrorCode ierr,(*f)(PC);

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

    528: /*@
    529:    PCILUReorderForNonzeroDiagonal - reorders rows/columns of matrix to remove zeros from diagonal

    531:    Collective on PC
    532:
    533:    Input Parameters:
    534: +  pc - the preconditioner context
    535: -  tol - diagonal entries smaller than this in absolute value are considered zero

    537:    Options Database Key:
    538: .  -pc_lu_nonzeros_along_diagonal

    540:    Level: intermediate

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

    544: .seealso: PCILUSetFill(), MatReorderForNonzeroDiagonal()
    545: @*/
    546: PetscErrorCode PETSCKSP_DLLEXPORT PCILUReorderForNonzeroDiagonal(PC pc,PetscReal rtol)
    547: {
    548:   PetscErrorCode ierr,(*f)(PC,PetscReal);

    552:   PetscObjectQueryFunction((PetscObject)pc,"PCILUReorderForNonzeroDiagonal_C",(void (**)(void))&f);
    553:   if (f) {
    554:     (*f)(pc,rtol);
    555:   }
    556:   return(0);
    557: }

    561: /*@
    562:     PCILUSetPivotInBlocks - Determines if pivoting is done while factoring each block
    563:       with BAIJ or SBAIJ matrices

    565:     Collective on PC

    567:     Input Parameters:
    568: +   pc - the preconditioner context
    569: -   pivot - PETSC_TRUE or PETSC_FALSE

    571:     Options Database Key:
    572: .   -pc_ilu_pivot_in_blocks <true,false>

    574:     Level: intermediate

    576: .seealso: PCIILUSetMatOrdering(), PCILUSetPivoting()
    577: @*/
    578: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetPivotInBlocks(PC pc,PetscTruth pivot)
    579: {
    580:   PetscErrorCode ierr,(*f)(PC,PetscTruth);

    583:   PetscObjectQueryFunction((PetscObject)pc,"PCILUSetPivotInBlocks_C",(void (**)(void))&f);
    584:   if (f) {
    585:     (*f)(pc,pivot);
    586:   }
    587:   return(0);
    588: }

    590: /* ------------------------------------------------------------------------------------------*/

    595: PetscErrorCode PETSCKSP_DLLEXPORT PCILUSetPivotInBlocks_ILU(PC pc,PetscTruth pivot)
    596: {
    597:   PC_ILU *dir = (PC_ILU*)pc->data;

    600:   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
    601:   return(0);
    602: }

    607: static PetscErrorCode PCSetFromOptions_ILU(PC pc)
    608: {
    610:   PetscInt       dtmax = 3,itmp;
    611:   PetscTruth     flg,set;
    612:   PetscReal      dt[3];
    613:   char           tname[256];
    614:   PC_ILU         *ilu = (PC_ILU*)pc->data;
    615:   PetscFList     ordlist;
    616:   PetscReal      tol;

    619:   MatOrderingRegisterAll(PETSC_NULL);
    620:   PetscOptionsHead("ILU Options");
    621:     PetscOptionsInt("-pc_ilu_levels","levels of fill","PCILUSetLevels",(PetscInt)ilu->info.levels,&itmp,&flg);
    622:     if (flg) ilu->info.levels = itmp;
    623:     PetscOptionsName("-pc_ilu_in_place","do factorization in place","PCILUSetUseInPlace",&ilu->inplace);
    624:     PetscOptionsName("-pc_ilu_diagonal_fill","Allow fill into empty diagonal entry","PCILUSetAllowDiagonalFill",&flg);
    625:     ilu->info.diagonal_fill = (double) flg;
    626:     PetscOptionsName("-pc_iludt_reuse_fill","Reuse fill from previous ILUdt","PCILUDTSetReuseFill",&ilu->reusefill);
    627:     PetscOptionsName("-pc_ilu_reuse_ordering","Reuse previous reordering","PCILUSetReuseOrdering",&ilu->reuseordering);

    629:     PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
    630:     if (flg) {
    631:       PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);
    632:     }
    633:     PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",ilu->info.shiftnz,&ilu->info.shiftnz,0);
    634:
    635:     PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",&flg);
    636:     if (flg) {
    637:       PetscOptionsInt("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",(PetscInt)ilu->info.shiftpd,&itmp,&flg);
    638:       if (flg && !itmp) {
    639:         PCFactorSetShiftPd(pc,PETSC_FALSE);
    640:       } else {
    641:         PCFactorSetShiftPd(pc,PETSC_TRUE);
    642:       }
    643:     }
    644:     PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",ilu->info.zeropivot,&ilu->info.zeropivot,0);

    646:     dt[0] = ilu->info.dt;
    647:     dt[1] = ilu->info.dtcol;
    648:     dt[2] = ilu->info.dtcount;
    649:     PetscOptionsRealArray("-pc_ilu_use_drop_tolerance","<dt,dtcol,maxrowcount>","PCILUSetUseDropTolerance",dt,&dtmax,&flg);
    650:     if (flg) {
    651:       PCILUSetUseDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
    652:     }
    653:     PetscOptionsReal("-pc_ilu_fill","Expected fill in factorization","PCILUSetFill",ilu->info.fill,&ilu->info.fill,&flg);
    654:     PetscOptionsName("-pc_ilu_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCILUReorderForNonzeroDiagonal",&flg);
    655:     if (flg) {
    656:       tol = PETSC_DECIDE;
    657:       PetscOptionsReal("-pc_ilu_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCILUReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);
    658:       PCILUReorderForNonzeroDiagonal(pc,tol);
    659:     }

    661:     MatGetOrderingList(&ordlist);
    662:     PetscOptionsList("-pc_ilu_mat_ordering_type","Reorder to reduce nonzeros in ILU","PCILUSetMatOrdering",ordlist,ilu->ordering,tname,256,&flg);
    663:     if (flg) {
    664:       PCILUSetMatOrdering(pc,tname);
    665:     }
    666:     flg = ilu->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
    667:     PetscOptionsTruth("-pc_ilu_pivot_in_blocks","Pivot inside matrix blocks for BAIJ and SBAIJ","PCILUSetPivotInBlocks",flg,&flg,&set);
    668:     if (set) {
    669:       PCILUSetPivotInBlocks(pc,flg);
    670:     }
    671:   PetscOptionsTail();
    672:   return(0);
    673: }

    677: static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
    678: {
    679:   PC_ILU         *ilu = (PC_ILU*)pc->data;
    681:   PetscTruth     isstring,iascii;

    684:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
    685:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    686:   if (iascii) {
    687:     if (ilu->usedt) {
    688:         PetscViewerASCIIPrintf(viewer,"  ILU: drop tolerance %g\n",ilu->info.dt);
    689:         PetscViewerASCIIPrintf(viewer,"  ILU: max nonzeros per row %D\n",(PetscInt)ilu->info.dtcount);
    690:         PetscViewerASCIIPrintf(viewer,"  ILU: column permutation tolerance %g\n",ilu->info.dtcol);
    691:     } else if (ilu->info.levels == 1) {
    692:         PetscViewerASCIIPrintf(viewer,"  ILU: %D level of fill\n",(PetscInt)ilu->info.levels);
    693:     } else {
    694:         PetscViewerASCIIPrintf(viewer,"  ILU: %D levels of fill\n",(PetscInt)ilu->info.levels);
    695:     }
    696:     PetscViewerASCIIPrintf(viewer,"  ILU: max fill ratio allocated %g\n",ilu->info.fill);
    697:     PetscViewerASCIIPrintf(viewer,"  ILU: tolerance for zero pivot %g\n",ilu->info.zeropivot);
    698:     if (ilu->info.shiftpd) {PetscViewerASCIIPrintf(viewer,"  ILU: using Manteuffel shift\n");}
    699:     if (ilu->inplace) {PetscViewerASCIIPrintf(viewer,"       in-place factorization\n");}
    700:     else              {PetscViewerASCIIPrintf(viewer,"       out-of-place factorization\n");}
    701:     PetscViewerASCIIPrintf(viewer,"       matrix ordering: %s\n",ilu->ordering);
    702:     if (ilu->reusefill)     {PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorization\n");}
    703:     if (ilu->reuseordering) {PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorization\n");}
    704:     if (ilu->fact) {
    705:       PetscViewerASCIIPrintf(viewer,"       Factored matrix follows\n");
    706:       PetscViewerASCIIPushTab(viewer);
    707:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
    708:       MatView(ilu->fact,viewer);
    709:       PetscViewerPopFormat(viewer);
    710:       PetscViewerASCIIPopTab(viewer);
    711:     }
    712:   } else if (isstring) {
    713:     PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)ilu->info.levels,ilu->ordering);
    714:   } else {
    715:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCILU",((PetscObject)viewer)->type_name);
    716:   }
    717:   return(0);
    718: }

    722: static PetscErrorCode PCSetUp_ILU(PC pc)
    723: {
    725:   PC_ILU         *ilu = (PC_ILU*)pc->data;

    728:   if (ilu->inplace) {
    729:     if (!pc->setupcalled) {

    731:       /* In-place factorization only makes sense with the natural ordering,
    732:          so we only need to get the ordering once, even if nonzero structure changes */
    733:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
    734:       if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
    735:       if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
    736:     }

    738:     /* In place ILU only makes sense with fill factor of 1.0 because
    739:        cannot have levels of fill */
    740:     ilu->info.fill          = 1.0;
    741:     ilu->info.diagonal_fill = 0;
    742:     MatILUFactor(pc->pmat,ilu->row,ilu->col,&ilu->info);
    743:     ilu->fact = pc->pmat;
    744:   } else if (ilu->usedt) {
    745:     if (!pc->setupcalled) {
    746:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
    747:       if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
    748:       if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
    749:       MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
    750:       PetscLogObjectParent(pc,ilu->fact);
    751:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
    752:       MatDestroy(ilu->fact);
    753:       if (!ilu->reuseordering) {
    754:         if (ilu->row) {ISDestroy(ilu->row);}
    755:         if (ilu->col) {ISDestroy(ilu->col);}
    756:         MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
    757:         if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
    758:         if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
    759:       }
    760:       MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
    761:       PetscLogObjectParent(pc,ilu->fact);
    762:     } else if (!ilu->reusefill) {
    763:       MatDestroy(ilu->fact);
    764:       MatILUDTFactor(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
    765:       PetscLogObjectParent(pc,ilu->fact);
    766:     } else {
    767:       MatLUFactorNumeric(pc->pmat,&ilu->info,&ilu->fact);
    768:     }
    769:    } else {
    770:     if (!pc->setupcalled) {
    771:       /* first time in so compute reordering and symbolic factorization */
    772:       MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
    773:       if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
    774:       if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
    775:       /*  Remove zeros along diagonal?     */
    776:       if (ilu->nonzerosalongdiagonal) {
    777:         MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
    778:       }
    779:       MatILUFactorSymbolic(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
    780:       PetscLogObjectParent(pc,ilu->fact);
    781:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
    782:       if (!ilu->reuseordering) {
    783:         /* compute a new ordering for the ILU */
    784:         ISDestroy(ilu->row);
    785:         ISDestroy(ilu->col);
    786:         MatGetOrdering(pc->pmat,ilu->ordering,&ilu->row,&ilu->col);
    787:         if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
    788:         if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
    789:         /*  Remove zeros along diagonal?     */
    790:         if (ilu->nonzerosalongdiagonal) {
    791:           MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
    792:         }
    793:       }
    794:       MatDestroy(ilu->fact);
    795:       MatILUFactorSymbolic(pc->pmat,ilu->row,ilu->col,&ilu->info,&ilu->fact);
    796:       PetscLogObjectParent(pc,ilu->fact);
    797:     }
    798:     MatLUFactorNumeric(pc->pmat,&ilu->info,&ilu->fact);
    799:   }
    800:   return(0);
    801: }

    805: static PetscErrorCode PCDestroy_ILU(PC pc)
    806: {
    807:   PC_ILU         *ilu = (PC_ILU*)pc->data;

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

    819: static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
    820: {
    821:   PC_ILU         *ilu = (PC_ILU*)pc->data;

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

    831: static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
    832: {
    833:   PC_ILU         *ilu = (PC_ILU*)pc->data;

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

    843: static PetscErrorCode PCGetFactoredMatrix_ILU(PC pc,Mat *mat)
    844: {
    845:   PC_ILU *ilu = (PC_ILU*)pc->data;

    848:   if (!ilu->fact) SETERRQ(PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
    849:   *mat = ilu->fact;
    850:   return(0);
    851: }

    853: /*MC
    854:      PCILU - Incomplete factorization preconditioners.

    856:    Options Database Keys:
    857: +  -pc_ilu_levels <k> - number of levels of fill for ILU(k)
    858: .  -pc_ilu_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
    859:                       its factorization (overwrites original matrix)
    860: .  -pc_ilu_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
    861: .  -pc_ilu_reuse_ordering - reuse ordering of factorized matrix from previous factorization
    862: .  -pc_ilu_use_drop_tolerance <dt,dtcol,maxrowcount> - use Saad's drop tolerance ILUdt
    863: .  -pc_ilu_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
    864: .  -pc_ilu_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
    865:                                    this decreases the chance of getting a zero pivot
    866: .  -pc_ilu_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
    867: .  -pc_ilu_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
    868:                              than 1 the diagonal blocks are factored with partial pivoting (this increases the
    869:                              stability of the ILU factorization
    870: .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
    871: -  -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
    872:    is optional with PETSC_TRUE being the default

    874:    Level: beginner

    876:   Concepts: incomplete factorization

    878:    Notes: Only implemented for some matrix formats. Not implemented in parallel

    880:           For BAIJ matrices this implements a point block ILU

    882: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
    883:            PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), PCILUSetUseDropTolerance(),
    884:            PCILUSetFill(), PCILUSetMatOrdering(), PCILUSetReuseOrdering(), PCILUDTSetReuseFill(),
    885:            PCILUSetLevels(), PCILUSetUseInPlace(), PCILUSetAllowDiagonalFill(), PCILUSetPivotInBlocks(),
    886:            PCFactorSetShiftNonzero(),PCFactorSetShiftPd()

    888: M*/

    893: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ILU(PC pc)
    894: {
    896:   PC_ILU         *ilu;

    899:   PetscNew(PC_ILU,&ilu);
    900:   PetscLogObjectMemory(pc,sizeof(PC_ILU));

    902:   ilu->fact                    = 0;
    903:   MatFactorInfoInitialize(&ilu->info);
    904:   ilu->info.levels             = 0;
    905:   ilu->info.fill               = 1.0;
    906:   ilu->col                     = 0;
    907:   ilu->row                     = 0;
    908:   ilu->inplace                 = PETSC_FALSE;
    909:   PetscStrallocpy(MATORDERING_NATURAL,&ilu->ordering);
    910:   ilu->reuseordering           = PETSC_FALSE;
    911:   ilu->usedt                   = PETSC_FALSE;
    912:   ilu->info.dt                 = PETSC_DEFAULT;
    913:   ilu->info.dtcount            = PETSC_DEFAULT;
    914:   ilu->info.dtcol              = PETSC_DEFAULT;
    915:   ilu->info.shiftnz            = 0.0;
    916:   ilu->info.shiftpd            = PETSC_FALSE;
    917:   ilu->info.shift_fraction     = 0.0;
    918:   ilu->info.zeropivot          = 1.e-12;
    919:   ilu->info.pivotinblocks      = 1.0;
    920:   ilu->reusefill               = PETSC_FALSE;
    921:   ilu->info.diagonal_fill      = 0;
    922:   pc->data                     = (void*)ilu;

    924:   pc->ops->destroy             = PCDestroy_ILU;
    925:   pc->ops->apply               = PCApply_ILU;
    926:   pc->ops->applytranspose      = PCApplyTranspose_ILU;
    927:   pc->ops->setup               = PCSetUp_ILU;
    928:   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
    929:   pc->ops->getfactoredmatrix   = PCGetFactoredMatrix_ILU;
    930:   pc->ops->view                = PCView_ILU;
    931:   pc->ops->applyrichardson     = 0;

    933:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_ILU",
    934:                     PCFactorSetZeroPivot_ILU);
    935:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_ILU",
    936:                     PCFactorSetShiftNonzero_ILU);
    937:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_ILU",
    938:                     PCFactorSetShiftPd_ILU);

    940:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetUseDropTolerance_C","PCILUSetUseDropTolerance_ILU",
    941:                     PCILUSetUseDropTolerance_ILU);
    942:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetFill_C","PCILUSetFill_ILU",
    943:                     PCILUSetFill_ILU);
    944:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetMatOrdering_C","PCILUSetMatOrdering_ILU",
    945:                     PCILUSetMatOrdering_ILU);
    946:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetReuseOrdering_C","PCILUSetReuseOrdering_ILU",
    947:                     PCILUSetReuseOrdering_ILU);
    948:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUDTSetReuseFill_C","PCILUDTSetReuseFill_ILUDT",
    949:                     PCILUDTSetReuseFill_ILUDT);
    950:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetLevels_C","PCILUSetLevels_ILU",
    951:                     PCILUSetLevels_ILU);
    952:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetUseInPlace_C","PCILUSetUseInPlace_ILU",
    953:                     PCILUSetUseInPlace_ILU);
    954:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetAllowDiagonalFill_C","PCILUSetAllowDiagonalFill_ILU",
    955:                     PCILUSetAllowDiagonalFill_ILU);
    956:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUSetPivotInBlocks_C","PCILUSetPivotInBlocks_ILU",
    957:                     PCILUSetPivotInBlocks_ILU);
    958:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCILUReorderForNonzeroDiagonal_C","PCILUReorderForNonzeroDiagonal_ILU",
    959:                     PCILUReorderForNonzeroDiagonal_ILU);
    960:   return(0);
    961: }