Actual source code: lu.c

      1: #define PETSCKSP_DLL

      3: /*
      4:    Defines a direct factorization preconditioner for any Mat implementation
      5:    Note: this need not be consided a preconditioner since it supplies
      6:          a direct solver.
      7: */

     9:  #include src/ksp/pc/pcimpl.h
     10:  #include src/ksp/pc/impls/factor/lu/lu.h

     15: PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetZeroPivot_LU(PC pc,PetscReal z)
     16: {
     17:   PC_LU *lu;

     20:   lu                 = (PC_LU*)pc->data;
     21:   lu->info.zeropivot = z;
     22:   return(0);
     23: }

     29: PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftNonzero_LU(PC pc,PetscReal shift)
     30: {
     31:   PC_LU *dir;

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

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

     62: PetscErrorCode PETSCKSP_DLLEXPORT PCLUReorderForNonzeroDiagonal_LU(PC pc,PetscReal z)
     63: {
     64:   PC_LU *lu = (PC_LU*)pc->data;

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

     80: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetReuseOrdering_LU(PC pc,PetscTruth flag)
     81: {
     82:   PC_LU *lu;

     85:   lu                = (PC_LU*)pc->data;
     86:   lu->reuseordering = flag;
     87:   return(0);
     88: }

     94: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetReuseFill_LU(PC pc,PetscTruth flag)
     95: {
     96:   PC_LU *lu;

     99:   lu = (PC_LU*)pc->data;
    100:   lu->reusefill = flag;
    101:   return(0);
    102: }

    107: static PetscErrorCode PCSetFromOptions_LU(PC pc)
    108: {
    109:   PC_LU          *lu = (PC_LU*)pc->data;
    111:   PetscTruth     flg,set;
    112:   char           tname[256];
    113:   PetscFList     ordlist;
    114:   PetscReal      tol;

    117:   MatOrderingRegisterAll(PETSC_NULL);
    118:   PetscOptionsHead("LU options");
    119:     PetscOptionsName("-pc_lu_in_place","Form LU in the same memory as the matrix","PCLUSetUseInPlace",&flg);
    120:     if (flg) {
    121:       PCLUSetUseInPlace(pc);
    122:     }
    123:     PetscOptionsReal("-pc_lu_fill","Expected non-zeros in LU/non-zeros in matrix","PCLUSetFill",lu->info.fill,&lu->info.fill,0);

    125:     PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
    126:     if (flg) {
    127:         PCFactorSetShiftNonzero(pc,(PetscReal) PETSC_DECIDE);
    128:     }
    129:     PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",lu->info.shiftnz,&lu->info.shiftnz,0);
    130:     PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",&flg);
    131:     if (flg) {
    132:       PCFactorSetShiftPd(pc,PETSC_TRUE);
    133:     }
    134:     PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",lu->info.zeropivot,&lu->info.zeropivot,0);

    136:     PetscOptionsName("-pc_lu_reuse_fill","Use fill from previous factorization","PCLUSetReuseFill",&flg);
    137:     if (flg) {
    138:       PCLUSetReuseFill(pc,PETSC_TRUE);
    139:     }
    140:     PetscOptionsName("-pc_lu_reuse_ordering","Reuse ordering from previous factorization","PCLUSetReuseOrdering",&flg);
    141:     if (flg) {
    142:       PCLUSetReuseOrdering(pc,PETSC_TRUE);
    143:     }

    145:     MatGetOrderingList(&ordlist);
    146:     PetscOptionsList("-pc_lu_mat_ordering_type","Reordering to reduce nonzeros in LU","PCLUSetMatOrdering",ordlist,lu->ordering,tname,256,&flg);
    147:     if (flg) {
    148:       PCLUSetMatOrdering(pc,tname);
    149:     }

    151:     PetscOptionsName("-pc_lu_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCLUReorderForNonzeroDiagonal",&flg);
    152:     if (flg) {
    153:       tol = PETSC_DECIDE;
    154:       PetscOptionsReal("-pc_lu_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCLUReorderForNonzeroDiagonal",lu->nonzerosalongdiagonaltol,&tol,0);
    155:       PCLUReorderForNonzeroDiagonal(pc,tol);
    156:     }

    158:     PetscOptionsReal("-pc_lu_pivoting","Pivoting tolerance (used only for some factorization)","PCLUSetPivoting",lu->info.dtcol,&lu->info.dtcol,&flg);

    160:     flg = lu->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
    161:     PetscOptionsTruth("-pc_lu_pivot_in_blocks","Pivot inside matrix blocks for BAIJ and SBAIJ","PCLUSetPivotInBlocks",flg,&flg,&set);
    162:     if (set) {
    163:       PCLUSetPivotInBlocks(pc,flg);
    164:     }
    165:   PetscOptionsTail();
    166:   return(0);
    167: }

    171: static PetscErrorCode PCView_LU(PC pc,PetscViewer viewer)
    172: {
    173:   PC_LU          *lu = (PC_LU*)pc->data;
    175:   PetscTruth     iascii,isstring;

    178:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    179:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
    180:   if (iascii) {
    181:     MatInfo info;

    183:     if (lu->inplace) {PetscViewerASCIIPrintf(viewer,"  LU: in-place factorization\n");}
    184:     else             {PetscViewerASCIIPrintf(viewer,"  LU: out-of-place factorization\n");}
    185:     PetscViewerASCIIPrintf(viewer,"    matrix ordering: %s\n",lu->ordering);
    186:     PetscViewerASCIIPrintf(viewer,"  LU: tolerance for zero pivot %g\n",lu->info.zeropivot);
    187:     if (lu->info.shiftpd) {PetscViewerASCIIPrintf(viewer,"  LU: using Manteuffel shift\n");}
    188:     if (lu->fact) {
    189:       MatGetInfo(lu->fact,MAT_LOCAL,&info);
    190:       PetscViewerASCIIPrintf(viewer,"    LU nonzeros %g\n",info.nz_used);
    191:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_FACTOR_INFO);
    192:       MatView(lu->fact,viewer);
    193:       PetscViewerPopFormat(viewer);
    194:     }
    195:     if (lu->reusefill)    {PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorization\n");}
    196:     if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorization\n");}
    197:   } else if (isstring) {
    198:     PetscViewerStringSPrintf(viewer," order=%s",lu->ordering);
    199:   } else {
    200:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCLU",((PetscObject)viewer)->type_name);
    201:   }
    202:   return(0);
    203: }

    207: static PetscErrorCode PCGetFactoredMatrix_LU(PC pc,Mat *mat)
    208: {
    209:   PC_LU *dir = (PC_LU*)pc->data;

    212:   if (!dir->fact) SETERRQ(PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
    213:   *mat = dir->fact;
    214:   return(0);
    215: }

    219: static PetscErrorCode PCSetUp_LU(PC pc)
    220: {
    222:   PC_LU          *dir = (PC_LU*)pc->data;

    225:   if (dir->reusefill && pc->setupcalled) dir->info.fill = dir->actualfill;

    227:   if (dir->inplace) {
    228:     if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(dir->row);}
    229:     if (dir->col) {ISDestroy(dir->col);}
    230:     MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
    231:     if (dir->row) {PetscLogObjectParent(pc,dir->row);  PetscLogObjectParent(pc,dir->col);}
    232:     MatLUFactor(pc->pmat,dir->row,dir->col,&dir->info);
    233:     dir->fact = pc->pmat;
    234:   } else {
    235:     MatInfo info;
    236:     if (!pc->setupcalled) {
    237:       MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
    238:       if (dir->nonzerosalongdiagonal) {
    239:         MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
    240:       }
    241:       if (dir->row) {PetscLogObjectParent(pc,dir->row); PetscLogObjectParent(pc,dir->col);}
    242:       MatLUFactorSymbolic(pc->pmat,dir->row,dir->col,&dir->info,&dir->fact);
    243:       MatGetInfo(dir->fact,MAT_LOCAL,&info);
    244:       dir->actualfill = info.fill_ratio_needed;
    245:       PetscLogObjectParent(pc,dir->fact);
    246:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
    247:       if (!dir->reuseordering) {
    248:         if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(dir->row);}
    249:         if (dir->col) {ISDestroy(dir->col);}
    250:         MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
    251:         if (dir->nonzerosalongdiagonal) {
    252:          MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
    253:         }
    254:         if (dir->row) {PetscLogObjectParent(pc,dir->row); PetscLogObjectParent(pc,dir->col);}
    255:       }
    256:       MatDestroy(dir->fact);
    257:       MatLUFactorSymbolic(pc->pmat,dir->row,dir->col,&dir->info,&dir->fact);
    258:       MatGetInfo(dir->fact,MAT_LOCAL,&info);
    259:       dir->actualfill = info.fill_ratio_needed;
    260:       PetscLogObjectParent(pc,dir->fact);
    261:     }
    262:     MatLUFactorNumeric(pc->pmat,&dir->info,&dir->fact);
    263:   }
    264:   return(0);
    265: }

    269: static PetscErrorCode PCDestroy_LU(PC pc)
    270: {
    271:   PC_LU          *dir = (PC_LU*)pc->data;

    275:   if (!dir->inplace && dir->fact) {MatDestroy(dir->fact);}
    276:   if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(dir->row);}
    277:   if (dir->col) {ISDestroy(dir->col);}
    278:   PetscStrfree(dir->ordering);
    279:   PetscFree(dir);
    280:   return(0);
    281: }

    285: static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
    286: {
    287:   PC_LU          *dir = (PC_LU*)pc->data;

    291:   if (dir->inplace) {MatSolve(pc->pmat,x,y);}
    292:   else              {MatSolve(dir->fact,x,y);}
    293:   return(0);
    294: }

    298: static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
    299: {
    300:   PC_LU          *dir = (PC_LU*)pc->data;

    304:   if (dir->inplace) {MatSolveTranspose(pc->pmat,x,y);}
    305:   else              {MatSolveTranspose(dir->fact,x,y);}
    306:   return(0);
    307: }

    309: /* -----------------------------------------------------------------------------------*/

    314: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetFill_LU(PC pc,PetscReal fill)
    315: {
    316:   PC_LU *dir;

    319:   dir = (PC_LU*)pc->data;
    320:   dir->info.fill = fill;
    321:   return(0);
    322: }

    328: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetUseInPlace_LU(PC pc)
    329: {
    330:   PC_LU *dir;

    333:   dir = (PC_LU*)pc->data;
    334:   dir->inplace = PETSC_TRUE;
    335:   return(0);
    336: }

    342: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetMatOrdering_LU(PC pc,MatOrderingType ordering)
    343: {
    344:   PC_LU          *dir = (PC_LU*)pc->data;

    348:   PetscStrfree(dir->ordering);
    349:   PetscStrallocpy(ordering,&dir->ordering);
    350:   return(0);
    351: }

    357: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetPivoting_LU(PC pc,PetscReal dtcol)
    358: {
    359:   PC_LU *dir = (PC_LU*)pc->data;

    362:   if (dtcol < 0.0 || dtcol > 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Column pivot tolerance is %g must be between 0 and 1",dtcol);
    363:   dir->info.dtcol = dtcol;
    364:   return(0);
    365: }

    371: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetPivotInBlocks_LU(PC pc,PetscTruth pivot)
    372: {
    373:   PC_LU *dir = (PC_LU*)pc->data;

    376:   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
    377:   return(0);
    378: }

    381: /* -----------------------------------------------------------------------------------*/

    385: /*@
    386:    PCLUReorderForNonzeroDiagonal - reorders rows/columns of matrix to remove zeros from diagonal

    388:    Collective on PC
    389:
    390:    Input Parameters:
    391: +  pc - the preconditioner context
    392: -  tol - diagonal entries smaller than this in absolute value are considered zero

    394:    Options Database Key:
    395: .  -pc_lu_nonzeros_along_diagonal

    397:    Level: intermediate

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

    401: .seealso: PCLUSetFill(), PCFactorSetShiftNonzero(), PCFactorSetZeroPivot(), MatReorderForNonzeroDiagonal()
    402: @*/
    403: PetscErrorCode PETSCKSP_DLLEXPORT PCLUReorderForNonzeroDiagonal(PC pc,PetscReal rtol)
    404: {
    405:   PetscErrorCode ierr,(*f)(PC,PetscReal);

    409:   PetscObjectQueryFunction((PetscObject)pc,"PCLUReorderForNonzeroDiagonal_C",(void (**)(void))&f);
    410:   if (f) {
    411:     (*f)(pc,rtol);
    412:   }
    413:   return(0);
    414: }

    418: /*@
    419:    PCLUSetReuseOrdering - When similar matrices are factored, this
    420:    causes the ordering computed in the first factor to be used for all
    421:    following factors; applies to both fill and drop tolerance LUs.

    423:    Collective on PC

    425:    Input Parameters:
    426: +  pc - the preconditioner context
    427: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

    429:    Options Database Key:
    430: .  -pc_lu_reuse_ordering - Activate PCLUSetReuseOrdering()

    432:    Level: intermediate

    434: .keywords: PC, levels, reordering, factorization, incomplete, LU

    436: .seealso: PCLUSetReuseFill(), PCILUSetReuseOrdering(), PCILUDTSetReuseFill()
    437: @*/
    438: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetReuseOrdering(PC pc,PetscTruth flag)
    439: {
    440:   PetscErrorCode ierr,(*f)(PC,PetscTruth);

    444:   PetscObjectQueryFunction((PetscObject)pc,"PCLUSetReuseOrdering_C",(void (**)(void))&f);
    445:   if (f) {
    446:     (*f)(pc,flag);
    447:   }
    448:   return(0);
    449: }

    453: /*@
    454:    PCLUSetReuseFill - When matrices with same nonzero structure are LU factored,
    455:    this causes later ones to use the fill computed in the initial factorization.

    457:    Collective on PC

    459:    Input Parameters:
    460: +  pc - the preconditioner context
    461: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

    463:    Options Database Key:
    464: .  -pc_lu_reuse_fill - Activates PCLUSetReuseFill()

    466:    Level: intermediate

    468: .keywords: PC, levels, reordering, factorization, incomplete, LU

    470: .seealso: PCILUSetReuseOrdering(), PCLUSetReuseOrdering(), PCILUDTSetReuseFill()
    471: @*/
    472: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetReuseFill(PC pc,PetscTruth flag)
    473: {
    474:   PetscErrorCode ierr,(*f)(PC,PetscTruth);

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

    487: /*@
    488:    PCLUSetFill - Indicate the amount of fill you expect in the factored matrix,
    489:    fill = number nonzeros in factor/number nonzeros in original matrix.

    491:    Collective on PC
    492:
    493:    Input Parameters:
    494: +  pc - the preconditioner context
    495: -  fill - amount of expected fill

    497:    Options Database Key:
    498: .  -pc_lu_fill <fill> - Sets fill amount

    500:    Level: intermediate

    502:    Note:
    503:    For sparse matrix factorizations it is difficult to predict how much
    504:    fill to expect. By running with the option -log_info PETSc will print the
    505:    actual amount of fill used; allowing you to set the value accurately for
    506:    future runs. Default PETSc uses a value of 5.0

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

    510: .seealso: PCILUSetFill()
    511: @*/
    512: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetFill(PC pc,PetscReal fill)
    513: {
    514:   PetscErrorCode ierr,(*f)(PC,PetscReal);

    518:   if (fill < 1.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Fill factor cannot be less then 1.0");
    519:   PetscObjectQueryFunction((PetscObject)pc,"PCLUSetFill_C",(void (**)(void))&f);
    520:   if (f) {
    521:     (*f)(pc,fill);
    522:   }
    523:   return(0);
    524: }

    528: /*@
    529:    PCLUSetUseInPlace - Tells the system to do an in-place factorization.
    530:    For dense matrices, this enables the solution of much larger problems.
    531:    For sparse matrices the factorization cannot be done truly in-place
    532:    so this does not save memory during the factorization, but after the matrix
    533:    is factored, the original unfactored matrix is freed, thus recovering that
    534:    space.

    536:    Collective on PC

    538:    Input Parameters:
    539: .  pc - the preconditioner context

    541:    Options Database Key:
    542: .  -pc_lu_in_place - Activates in-place factorization

    544:    Notes:
    545:    PCLUSetUseInplace() can only be used with the KSP method KSPPREONLY or when
    546:    a different matrix is provided for the multiply and the preconditioner in
    547:    a call to KSPSetOperators().
    548:    This is because the Krylov space methods require an application of the
    549:    matrix multiplication, which is not possible here because the matrix has
    550:    been factored in-place, replacing the original matrix.

    552:    Level: intermediate

    554: .keywords: PC, set, factorization, direct, inplace, in-place, LU

    556: .seealso: PCILUSetUseInPlace()
    557: @*/
    558: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetUseInPlace(PC pc)
    559: {
    560:   PetscErrorCode ierr,(*f)(PC);

    564:   PetscObjectQueryFunction((PetscObject)pc,"PCLUSetUseInPlace_C",(void (**)(void))&f);
    565:   if (f) {
    566:     (*f)(pc);
    567:   }
    568:   return(0);
    569: }

    573: /*@C
    574:     PCLUSetMatOrdering - Sets the ordering routine (to reduce fill) to
    575:     be used in the LU factorization.

    577:     Collective on PC

    579:     Input Parameters:
    580: +   pc - the preconditioner context
    581: -   ordering - the matrix ordering name, for example, MATORDERING_ND or MATORDERING_RCM

    583:     Options Database Key:
    584: .   -pc_lu_mat_ordering_type <nd,rcm,...> - Sets ordering routine

    586:     Level: intermediate

    588:     Notes: nested dissection is used by default

    590: .seealso: PCILUSetMatOrdering()
    591: @*/
    592: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetMatOrdering(PC pc,MatOrderingType ordering)
    593: {
    594:   PetscErrorCode ierr,(*f)(PC,MatOrderingType);

    597:   PetscObjectQueryFunction((PetscObject)pc,"PCLUSetMatOrdering_C",(void (**)(void))&f);
    598:   if (f) {
    599:     (*f)(pc,ordering);
    600:   }
    601:   return(0);
    602: }

    606: /*@
    607:     PCLUSetPivoting - Determines when pivoting is done during LU.
    608:       For PETSc dense matrices column pivoting is always done, for PETSc sparse matrices
    609:       it is never done. For the Matlab and SuperLU factorization this is used.

    611:     Collective on PC

    613:     Input Parameters:
    614: +   pc - the preconditioner context
    615: -   dtcol - 0.0 implies no pivoting, 1.0 complete pivoting (slower, requires more memory but more stable)

    617:     Options Database Key:
    618: .   -pc_lu_pivoting <dtcol>

    620:     Level: intermediate

    622: .seealso: PCILUSetMatOrdering(), PCLUSetPivotInBlocks()
    623: @*/
    624: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetPivoting(PC pc,PetscReal dtcol)
    625: {
    626:   PetscErrorCode ierr,(*f)(PC,PetscReal);

    629:   PetscObjectQueryFunction((PetscObject)pc,"PCLUSetPivoting_C",(void (**)(void))&f);
    630:   if (f) {
    631:     (*f)(pc,dtcol);
    632:   }
    633:   return(0);
    634: }

    638: /*@
    639:     PCLUSetPivotInBlocks - Determines if pivoting is done while factoring each block
    640:       with BAIJ or SBAIJ matrices

    642:     Collective on PC

    644:     Input Parameters:
    645: +   pc - the preconditioner context
    646: -   pivot - PETSC_TRUE or PETSC_FALSE

    648:     Options Database Key:
    649: .   -pc_lu_pivot_in_blocks <true,false>

    651:     Level: intermediate

    653: .seealso: PCILUSetMatOrdering(), PCLUSetPivoting()
    654: @*/
    655: PetscErrorCode PETSCKSP_DLLEXPORT PCLUSetPivotInBlocks(PC pc,PetscTruth pivot)
    656: {
    657:   PetscErrorCode ierr,(*f)(PC,PetscTruth);

    660:   PetscObjectQueryFunction((PetscObject)pc,"PCLUSetPivotInBlocks_C",(void (**)(void))&f);
    661:   if (f) {
    662:     (*f)(pc,pivot);
    663:   }
    664:   return(0);
    665: }

    667: /* ------------------------------------------------------------------------ */

    669: /*MC
    670:    PCLU - Uses a direct solver, based on LU factorization, as a preconditioner

    672:    Options Database Keys:
    673: +  -pc_lu_reuse_ordering - Activate PCLUSetReuseOrdering()
    674: .  -pc_lu_reuse_fill - Activates PCLUSetReuseFill()
    675: .  -pc_lu_fill <fill> - Sets fill amount
    676: .  -pc_lu_in_place - Activates in-place factorization
    677: .  -pc_lu_mat_ordering_type <nd,rcm,...> - Sets ordering routine
    678: .  -pc_lu_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
    679:                                          stability of factorization.
    680: .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
    681: -  -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
    682:    is optional with PETSC_TRUE being the default

    684:    Notes: Not all options work for all matrix formats
    685:           Run with -help to see additional options for particular matrix formats or factorization
    686:           algorithms

    688:    Level: beginner

    690:    Concepts: LU factorization, direct solver

    692:    Notes: Usually this will compute an "exact" solution in one iteration and does
    693:           not need a Krylov method (i.e. you can use -ksp_type preonly, or
    694:           KSPSetType(ksp,KSPPREONLY) for the Krylov method

    696: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    697:            PCILU, PCCHOLESKY, PCICC, PCLUSetReuseOrdering(), PCLUSetReuseFill(), PCGetFactoredMatrix(),
    698:            PCLUSetFill(), PCLUSetUseInPlace(), PCLUSetMatOrdering(), PCFactorSetPivoting(),
    699:            PCLUSetPivotingInBlocks(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd()
    700: M*/

    705: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_LU(PC pc)
    706: {
    708:   PetscMPIInt    size;
    709:   PC_LU          *dir;

    712:   PetscNew(PC_LU,&dir);
    713:   PetscLogObjectMemory(pc,sizeof(PC_LU));

    715:   MatFactorInfoInitialize(&dir->info);
    716:   dir->fact                  = 0;
    717:   dir->inplace               = PETSC_FALSE;
    718:   dir->nonzerosalongdiagonal = PETSC_FALSE;

    720:   dir->info.fill           = 5.0;
    721:   dir->info.dtcol          = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
    722:   dir->info.shiftnz        = 0.0;
    723:   dir->info.zeropivot      = 1.e-12;
    724:   dir->info.pivotinblocks  = 1.0;
    725:   dir->info.shiftpd        = PETSC_FALSE;
    726:   dir->info.shift_fraction = 0.0;
    727:   dir->col                 = 0;
    728:   dir->row                 = 0;
    729:   MPI_Comm_size(pc->comm,&size);
    730:   if (size == 1) {
    731:     PetscStrallocpy(MATORDERING_ND,&dir->ordering);
    732:   } else {
    733:     PetscStrallocpy(MATORDERING_NATURAL,&dir->ordering);
    734:   }
    735:   dir->reusefill        = PETSC_FALSE;
    736:   dir->reuseordering    = PETSC_FALSE;
    737:   pc->data              = (void*)dir;

    739:   pc->ops->destroy           = PCDestroy_LU;
    740:   pc->ops->apply             = PCApply_LU;
    741:   pc->ops->applytranspose    = PCApplyTranspose_LU;
    742:   pc->ops->setup             = PCSetUp_LU;
    743:   pc->ops->setfromoptions    = PCSetFromOptions_LU;
    744:   pc->ops->view              = PCView_LU;
    745:   pc->ops->applyrichardson   = 0;
    746:   pc->ops->getfactoredmatrix = PCGetFactoredMatrix_LU;

    748:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_LU",
    749:                     PCFactorSetZeroPivot_LU);
    750:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_LU",
    751:                     PCFactorSetShiftNonzero_LU);
    752:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_LU",
    753:                     PCFactorSetShiftPd_LU);

    755:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCLUSetFill_C","PCLUSetFill_LU",
    756:                     PCLUSetFill_LU);
    757:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCLUSetUseInPlace_C","PCLUSetUseInPlace_LU",
    758:                     PCLUSetUseInPlace_LU);
    759:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCLUSetMatOrdering_C","PCLUSetMatOrdering_LU",
    760:                     PCLUSetMatOrdering_LU);
    761:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCLUSetReuseOrdering_C","PCLUSetReuseOrdering_LU",
    762:                     PCLUSetReuseOrdering_LU);
    763:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCLUSetReuseFill_C","PCLUSetReuseFill_LU",
    764:                     PCLUSetReuseFill_LU);
    765:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCLUSetPivoting_C","PCLUSetPivoting_LU",
    766:                     PCLUSetPivoting_LU);
    767:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCLUSetPivotInBlocks_C","PCLUSetPivotInBlocks_LU",
    768:                     PCLUSetPivotInBlocks_LU);
    769:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCLUReorderForNonzeroDiagonal_C","PCLUReorderForNonzeroDiagonal_LU",
    770:                     PCLUReorderForNonzeroDiagonal_LU);
    771:   return(0);
    772: }