Actual source code: icc.c

      2: #include <../src/ksp/pc/impls/factor/icc/icc.h>

      4: static PetscErrorCode PCSetUp_ICC(PC pc)
      5: {
      6:   PC_ICC                 *icc = (PC_ICC*)pc->data;
      7:   IS                     perm = NULL,cperm = NULL;
      8:   PetscErrorCode         ierr;
      9:   MatInfo                info;
     10:   MatSolverType          stype;
     11:   MatFactorError         err;

     14:   pc->failedreason = PC_NOERROR;

     16:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
     17:   if (!pc->setupcalled) {
     18:     if (!((PC_Factor*)icc)->fact) {
     19:       PetscBool useordering;
     20:       MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
     21:       MatFactorGetUseOrdering(((PC_Factor*)icc)->fact,&useordering);
     22:       if (useordering) {
     23:         MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
     24:       }
     25:     }
     26:     MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
     27:   } else if (pc->flag != SAME_NONZERO_PATTERN) {
     28:     PetscBool useordering;
     29:     MatDestroy(&((PC_Factor*)icc)->fact);
     30:     MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
     31:     MatFactorGetUseOrdering(((PC_Factor*)icc)->fact,&useordering);
     32:     if (useordering) {
     33:       MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
     34:     }
     35:     MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
     36:   }
     37:   MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
     38:   icc->hdr.actualfill = info.fill_ratio_needed;

     40:   ISDestroy(&cperm);
     41:   ISDestroy(&perm);

     43:   MatFactorGetError(((PC_Factor*)icc)->fact,&err);
     44:   if (err) { /* FactorSymbolic() fails */
     45:     pc->failedreason = (PCFailedReason)err;
     46:     return(0);
     47:   }

     49:   MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
     50:   MatFactorGetError(((PC_Factor*)icc)->fact,&err);
     51:   if (err) { /* FactorNumeric() fails */
     52:     pc->failedreason = (PCFailedReason)err;
     53:   }

     55:   PCFactorGetMatSolverType(pc,&stype);
     56:   if (!stype) {
     57:     MatSolverType solverpackage;
     58:     MatFactorGetSolverType(((PC_Factor*)icc)->fact,&solverpackage);
     59:     PCFactorSetMatSolverType(pc,solverpackage);
     60:   }
     61:   return(0);
     62: }

     64: static PetscErrorCode PCReset_ICC(PC pc)
     65: {
     66:   PC_ICC         *icc = (PC_ICC*)pc->data;

     70:   MatDestroy(&((PC_Factor*)icc)->fact);
     71:   return(0);
     72: }

     74: static PetscErrorCode PCDestroy_ICC(PC pc)
     75: {
     76:   PC_ICC         *icc = (PC_ICC*)pc->data;

     80:   PCReset_ICC(pc);
     81:   PetscFree(((PC_Factor*)icc)->ordering);
     82:   PetscFree(((PC_Factor*)icc)->solvertype);
     83:   PetscFree(pc->data);
     84:   return(0);
     85: }

     87: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
     88: {
     89:   PC_ICC         *icc = (PC_ICC*)pc->data;

     93:   MatSolve(((PC_Factor*)icc)->fact,x,y);
     94:   return(0);
     95: }

     97: static PetscErrorCode PCMatApply_ICC(PC pc,Mat X,Mat Y)
     98: {
     99:   PC_ICC         *icc = (PC_ICC*)pc->data;

    103:   MatMatSolve(((PC_Factor*)icc)->fact,X,Y);
    104:   return(0);
    105: }

    107: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
    108: {
    110:   PC_ICC         *icc = (PC_ICC*)pc->data;

    113:   MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
    114:   return(0);
    115: }

    117: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
    118: {
    120:   PC_ICC         *icc = (PC_ICC*)pc->data;

    123:   MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
    124:   return(0);
    125: }

    127: static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc)
    128: {
    129:   PC_ICC         *icc = (PC_ICC*)pc->data;
    130:   PetscBool      flg;
    132:   /* PetscReal      dt[3];*/

    135:   PetscOptionsHead(PetscOptionsObject,"ICC Options");
    136:   PCSetFromOptions_Factor(PetscOptionsObject,pc);

    138:   PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
    139:   /*dt[0] = ((PC_Factor*)icc)->info.dt;
    140:   dt[1] = ((PC_Factor*)icc)->info.dtcol;
    141:   dt[2] = ((PC_Factor*)icc)->info.dtcount;
    142:   PetscInt       dtmax = 3;
    143:   PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
    144:   if (flg) {
    145:     PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
    146:   }
    147:   */
    148:   PetscOptionsTail();
    149:   return(0);
    150: }

    152: extern PetscErrorCode  PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);

    154: /*MC
    155:      PCICC - Incomplete Cholesky factorization preconditioners.

    157:    Options Database Keys:
    158: +  -pc_factor_levels <k> - number of levels of fill for ICC(k)
    159: .  -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
    160:                       its factorization (overwrites original matrix)
    161: .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
    162: -  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix

    164:    Level: beginner

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

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

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

    173:           By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE);
    174:           to turn off the shift.

    176:    References:
    177: .  1. - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
    178:       Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
    179:       Science and Engineering, Kluwer.

    182: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
    183:            PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
    184:            PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
    185:            PCFactorSetLevels()

    187: M*/

    189: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
    190: {
    192:   PC_ICC         *icc;

    195:   PetscNewLog(pc,&icc);
    196:   pc->data = (void*)icc;
    197:   PCFactorInitialize(pc);
    198:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);

    200:   ((PC_Factor*)icc)->factortype     = MAT_FACTOR_ICC;
    201:   ((PC_Factor*)icc)->info.fill      = 1.0;
    202:   ((PC_Factor*)icc)->info.dtcol     = PETSC_DEFAULT;
    203:   ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;

    205:   pc->ops->apply               = PCApply_ICC;
    206:   pc->ops->matapply            = PCMatApply_ICC;
    207:   pc->ops->applytranspose      = PCApply_ICC;
    208:   pc->ops->setup               = PCSetUp_ICC;
    209:   pc->ops->reset               = PCReset_ICC;
    210:   pc->ops->destroy             = PCDestroy_ICC;
    211:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
    212:   pc->ops->view                = PCView_Factor;
    213:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
    214:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
    215:   return(0);
    216: }