Actual source code: icc.c

      1: #define PETSCKSP_DLL

      3: /*
      4:    Defines a Cholesky factorization preconditioner for any Mat implementation.
      5:   Presently only provided for MPIRowbs format (i.e. BlockSolve).
      6: */

     8:  #include src/ksp/pc/impls/factor/icc/icc.h

     13: PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetZeroPivot_ICC(PC pc,PetscReal z)
     14: {
     15:   PC_ICC *icc;

     18:   icc                 = (PC_ICC*)pc->data;
     19:   icc->info.zeropivot = z;
     20:   return(0);
     21: }

     27: PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftNonzero_ICC(PC pc,PetscReal shift)
     28: {
     29:   PC_ICC *dir;

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

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

     60: PetscErrorCode PETSCKSP_DLLEXPORT PCICCSetMatOrdering_ICC(PC pc,MatOrderingType ordering)
     61: {
     62:   PC_ICC         *dir = (PC_ICC*)pc->data;
     64:
     66:   PetscStrfree(dir->ordering);
     67:   PetscStrallocpy(ordering,&dir->ordering);
     68:   return(0);
     69: }

     75: PetscErrorCode PETSCKSP_DLLEXPORT PCICCSetFill_ICC(PC pc,PetscReal fill)
     76: {
     77:   PC_ICC *dir;

     80:   dir            = (PC_ICC*)pc->data;
     81:   dir->info.fill = fill;
     82:   return(0);
     83: }

     89: PetscErrorCode PETSCKSP_DLLEXPORT PCICCSetLevels_ICC(PC pc,PetscInt levels)
     90: {
     91:   PC_ICC *icc;

     94:   icc = (PC_ICC*)pc->data;
     95:   icc->info.levels = levels;
     96:   return(0);
     97: }

    102: /*@
    103:     PCICCSetMatOrdering - Sets the ordering routine (to reduce fill) to
    104:     be used it the ICC factorization.

    106:     Collective on PC

    108:     Input Parameters:
    109: +   pc - the preconditioner context
    110: -   ordering - the matrix ordering name, for example, MATORDERING_ND or MATORDERING_RCM

    112:     Options Database Key:
    113: .   -pc_icc_mat_ordering_type <nd,rcm,...> - Sets ordering routine

    115:     Level: intermediate

    117: .seealso: PCLUSetMatOrdering()

    119: .keywords: PC, ICC, set, matrix, reordering

    121: @*/
    122: PetscErrorCode PETSCKSP_DLLEXPORT PCICCSetMatOrdering(PC pc,MatOrderingType ordering)
    123: {
    124:   PetscErrorCode ierr,(*f)(PC,MatOrderingType);

    128:   PetscObjectQueryFunction((PetscObject)pc,"PCICCSetMatOrdering_C",(void (**)(void))&f);
    129:   if (f) {
    130:     (*f)(pc,ordering);
    131:   }
    132:   return(0);
    133: }

    137: /*@
    138:    PCICCSetLevels - Sets the number of levels of fill to use.

    140:    Collective on PC

    142:    Input Parameters:
    143: +  pc - the preconditioner context
    144: -  levels - number of levels of fill

    146:    Options Database Key:
    147: .  -pc_icc_levels <levels> - Sets fill level

    149:    Level: intermediate

    151:    Concepts: ICC^setting levels of fill

    153: @*/
    154: PetscErrorCode PETSCKSP_DLLEXPORT PCICCSetLevels(PC pc,PetscInt levels)
    155: {
    156:   PetscErrorCode ierr,(*f)(PC,PetscInt);

    160:   if (levels < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"negative levels");
    161:   PetscObjectQueryFunction((PetscObject)pc,"PCICCSetLevels_C",(void (**)(void))&f);
    162:   if (f) {
    163:     (*f)(pc,levels);
    164:   }
    165:   return(0);
    166: }

    170: /*@
    171:    PCICCSetFill - Indicate the amount of fill you expect in the factored matrix,
    172:    where fill = number nonzeros in factor/number nonzeros in original matrix.

    174:    Collective on PC

    176:    Input Parameters:
    177: +  pc - the preconditioner context
    178: -  fill - amount of expected fill

    180:    Options Database Key:
    181: $  -pc_icc_fill <fill>

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

    189:    Level: intermediate

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

    193: .seealso: PCLUSetFill()
    194: @*/
    195: PetscErrorCode PETSCKSP_DLLEXPORT PCICCSetFill(PC pc,PetscReal fill)
    196: {
    197:   PetscErrorCode ierr,(*f)(PC,PetscReal);

    201:   if (fill < 1.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Fill factor cannot be less than 1.0");
    202:   PetscObjectQueryFunction((PetscObject)pc,"PCICCSetFill_C",(void (**)(void))&f);
    203:   if (f) {
    204:     (*f)(pc,fill);
    205:   }
    206:   return(0);
    207: }

    211: static PetscErrorCode PCSetup_ICC(PC pc)
    212: {
    213:   PC_ICC         *icc = (PC_ICC*)pc->data;
    214:   IS             perm,cperm;

    218:   MatGetOrdering(pc->pmat,icc->ordering,&perm,&cperm);

    220:   if (!pc->setupcalled) {
    221:     MatICCFactorSymbolic(pc->pmat,perm,&icc->info,&icc->fact);
    222:   } else if (pc->flag != SAME_NONZERO_PATTERN) {
    223:     MatDestroy(icc->fact);
    224:     MatICCFactorSymbolic(pc->pmat,perm,&icc->info,&icc->fact);
    225:   }
    226:   ISDestroy(cperm);
    227:   ISDestroy(perm);
    228:   MatCholeskyFactorNumeric(pc->pmat,&icc->info,&icc->fact);
    229:   return(0);
    230: }

    234: static PetscErrorCode PCDestroy_ICC(PC pc)
    235: {
    236:   PC_ICC         *icc = (PC_ICC*)pc->data;

    240:   if (icc->fact) {MatDestroy(icc->fact);}
    241:   PetscStrfree(icc->ordering);
    242:   PetscFree(icc);
    243:   return(0);
    244: }

    248: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
    249: {
    250:   PC_ICC         *icc = (PC_ICC*)pc->data;

    254:   MatSolve(icc->fact,x,y);
    255:   return(0);
    256: }

    260: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
    261: {
    263:   PC_ICC         *icc = (PC_ICC*)pc->data;

    266:   MatForwardSolve(icc->fact,x,y);
    267:   return(0);
    268: }

    272: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
    273: {
    275:   PC_ICC         *icc = (PC_ICC*)pc->data;

    278:   MatBackwardSolve(icc->fact,x,y);
    279:   return(0);
    280: }

    284: static PetscErrorCode PCGetFactoredMatrix_ICC(PC pc,Mat *mat)
    285: {
    286:   PC_ICC *icc = (PC_ICC*)pc->data;

    289:   *mat = icc->fact;
    290:   return(0);
    291: }

    295: static PetscErrorCode PCSetFromOptions_ICC(PC pc)
    296: {
    297:   PC_ICC         *icc = (PC_ICC*)pc->data;
    298:   char           tname[256];
    299:   PetscTruth     flg;
    301:   PetscFList     ordlist;

    304:   MatOrderingRegisterAll(PETSC_NULL);
    305:   PetscOptionsHead("ICC Options");
    306:     PetscOptionsReal("-pc_icc_levels","levels of fill","PCICCSetLevels",icc->info.levels,&icc->info.levels,&flg);
    307:     PetscOptionsReal("-pc_icc_fill","Expected fill in factorization","PCICCSetFill",icc->info.fill,&icc->info.fill,&flg);
    308:     MatGetOrderingList(&ordlist);
    309:     PetscOptionsList("-pc_icc_mat_ordering_type","Reorder to reduce nonzeros in ICC","PCICCSetMatOrdering",ordlist,icc->ordering,tname,256,&flg);
    310:     if (flg) {
    311:       PCICCSetMatOrdering(pc,tname);
    312:     }
    313:     PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
    314:     if (flg) {
    315:       PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);
    316:     }
    317:     PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",icc->info.shiftnz,&icc->info.shiftnz,0);
    318:     PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCICCSetShift",&flg);
    319:     if (flg) {
    320:       PCFactorSetShiftPd(pc,PETSC_TRUE);
    321:     } else {
    322:       PCFactorSetShiftPd(pc,PETSC_FALSE);
    323:     }
    324:     PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",icc->info.zeropivot,&icc->info.zeropivot,0);
    325:
    326:   PetscOptionsTail();
    327:   return(0);
    328: }

    332: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
    333: {
    334:   PC_ICC         *icc = (PC_ICC*)pc->data;
    336:   PetscTruth     isstring,iascii;

    339:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
    340:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    341:   if (iascii) {
    342:     if (icc->info.levels == 1) {
    343:         PetscViewerASCIIPrintf(viewer,"  ICC: %D level of fill\n",(PetscInt)icc->info.levels);
    344:     } else {
    345:         PetscViewerASCIIPrintf(viewer,"  ICC: %D levels of fill\n",(PetscInt)icc->info.levels);
    346:     }
    347:     PetscViewerASCIIPrintf(viewer,"  ICC: max fill ratio allocated %g\n",icc->info.fill);
    348:     if (icc->info.shiftpd) {PetscViewerASCIIPrintf(viewer,"  ICC: using Manteuffel shift\n");}
    349:   } else if (isstring) {
    350:     PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)icc->info.levels);
    351:   } else {
    352:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name);
    353:   }
    354:   return(0);
    355: }

    357: /*MC
    358:      PCICC - Incomplete Cholesky factorization preconditioners.

    360:    Options Database Keys:
    361: +  -pc_icc_levels <k> - number of levels of fill for ICC(k)
    362: .  -pc_icc_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
    363:                       its factorization (overwrites original matrix)
    364: .  -pc_icc_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
    365: .  -pc_icc_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
    366: .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
    367: -  -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
    368:    is optional with PETSC_TRUE being the default

    370:    Level: beginner

    372:   Concepts: incomplete Cholesky factorization

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

    376:           For BAIJ matrices this implements a point block ICC.

    378:           The Manteuffel shift is only implemented for matrices with block size 1

    380:           By default, the Manteuffel is applied (for matrices with block size 1). Call PCICCSetShift(pc,PETSC_FALSE);
    381:           to turn off the shift.

    384: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
    385:            PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(),
    386:            PCICCSetFill(), PCICCSetMatOrdering(), PCICCSetReuseOrdering(),
    387:            PCICCSetLevels(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd(),

    389: M*/

    394: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ICC(PC pc)
    395: {
    397:   PC_ICC         *icc;

    400:   PetscNew(PC_ICC,&icc);
    401:   PetscLogObjectMemory(pc,sizeof(PC_ICC));

    403:   icc->fact                  = 0;
    404:   PetscStrallocpy(MATORDERING_NATURAL,&icc->ordering);
    405:   MatFactorInfoInitialize(&icc->info);
    406:   icc->info.levels          = 0;
    407:   icc->info.fill          = 1.0;
    408:   icc->implctx            = 0;

    410:   icc->info.dtcol              = PETSC_DEFAULT;
    411:   icc->info.shiftnz            = 0.0;
    412:   icc->info.shiftpd            = PETSC_TRUE;
    413:   icc->info.shift_fraction     = 0.0;
    414:   icc->info.zeropivot          = 1.e-12;
    415:   pc->data                       = (void*)icc;

    417:   pc->ops->apply               = PCApply_ICC;
    418:   pc->ops->setup               = PCSetup_ICC;
    419:   pc->ops->destroy               = PCDestroy_ICC;
    420:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
    421:   pc->ops->view                = PCView_ICC;
    422:   pc->ops->getfactoredmatrix   = PCGetFactoredMatrix_ICC;
    423:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
    424:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;

    426:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_ICC",
    427:                     PCFactorSetZeroPivot_ICC);
    428:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_ICC",
    429:                     PCFactorSetShiftNonzero_ICC);
    430:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_ICC",
    431:                     PCFactorSetShiftPd_ICC);

    433:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCICCSetLevels_C","PCICCSetLevels_ICC",
    434:                     PCICCSetLevels_ICC);
    435:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCICCSetFill_C","PCICCSetFill_ICC",
    436:                     PCICCSetFill_ICC);
    437:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCICCSetMatOrdering_C","PCICCSetMatOrdering_ICC",
    438:                     PCICCSetMatOrdering_ICC);
    439:   return(0);
    440: }