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:   MatInfo        info;
  9:   MatSolverType  stype;
 10:   MatFactorError err;
 11:   PetscBool      canuseordering;
 12:   const char    *prefix;

 14:   pc->failedreason = PC_NOERROR;

 16:   PCGetOptionsPrefix(pc, &prefix);
 17:   MatSetOptionsPrefixFactor(pc->pmat, prefix);

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

 42:   ISDestroy(&cperm);
 43:   ISDestroy(&perm);

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

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

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

 66: static PetscErrorCode PCReset_ICC(PC pc)
 67: {
 68:   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;

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

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

 90:   MatSolve(((PC_Factor *)icc)->fact, x, y);
 91:   return 0;
 92: }

 94: static PetscErrorCode PCMatApply_ICC(PC pc, Mat X, Mat Y)
 95: {
 96:   PC_ICC *icc = (PC_ICC *)pc->data;

 98:   MatMatSolve(((PC_Factor *)icc)->fact, X, Y);
 99:   return 0;
100: }

102: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc, Vec x, Vec y)
103: {
104:   PC_ICC *icc = (PC_ICC *)pc->data;

106:   MatForwardSolve(((PC_Factor *)icc)->fact, x, y);
107:   return 0;
108: }

110: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc, Vec x, Vec y)
111: {
112:   PC_ICC *icc = (PC_ICC *)pc->data;

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

118: static PetscErrorCode PCSetFromOptions_ICC(PC pc, PetscOptionItems *PetscOptionsObject)
119: {
120:   PC_ICC   *icc = (PC_ICC *)pc->data;
121:   PetscBool flg;
122:   /* PetscReal      dt[3];*/

124:   PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options");
125:   PCSetFromOptions_Factor(pc, PetscOptionsObject);

127:   PetscOptionsReal("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", ((PC_Factor *)icc)->info.levels, &((PC_Factor *)icc)->info.levels, &flg);
128:   /*dt[0] = ((PC_Factor*)icc)->info.dt;
129:   dt[1] = ((PC_Factor*)icc)->info.dtcol;
130:   dt[2] = ((PC_Factor*)icc)->info.dtcount;
131:   PetscInt       dtmax = 3;
132:   PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
133:   if (flg) {
134:     PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
135:   }
136:   */
137:   PetscOptionsHeadEnd();
138:   return 0;
139: }

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

143: /*MC
144:      PCICC - Incomplete Cholesky factorization preconditioners.

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

153:    Level: beginner

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

158:    For `MATSEQBAIJ` matrices this implements a point block ICC.

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

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

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

170: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCILU`, `PCLU`, `PCCHOLESKY`,
171:           `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
172:           `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
173:           `PCFactorSetLevels()`
174: M*/

176: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
177: {
178:   PC_ICC *icc;

180:   PetscNew(&icc);
181:   pc->data = (void *)icc;
182:   PCFactorInitialize(pc, MAT_FACTOR_ICC);

184:   ((PC_Factor *)icc)->info.fill      = 1.0;
185:   ((PC_Factor *)icc)->info.dtcol     = PETSC_DEFAULT;
186:   ((PC_Factor *)icc)->info.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;

188:   pc->ops->apply               = PCApply_ICC;
189:   pc->ops->matapply            = PCMatApply_ICC;
190:   pc->ops->applytranspose      = PCApply_ICC;
191:   pc->ops->setup               = PCSetUp_ICC;
192:   pc->ops->reset               = PCReset_ICC;
193:   pc->ops->destroy             = PCDestroy_ICC;
194:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
195:   pc->ops->view                = PCView_Factor;
196:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
197:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
198:   return 0;
199: }