Actual source code: icc.c

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

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

 13:   PetscFunctionBegin;
 14:   pc->failedreason = PC_NOERROR;

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

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

 41:   PetscCall(ISDestroy(&cperm));
 42:   PetscCall(ISDestroy(&perm));

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

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

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

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

 69:   PetscFunctionBegin;
 70:   PetscCall(MatDestroy(&((PC_Factor *)icc)->fact));
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

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

 78:   PetscFunctionBegin;
 79:   PetscCall(PCReset_ICC(pc));
 80:   PetscCall(PetscFree(((PC_Factor *)icc)->ordering));
 81:   PetscCall(PetscFree(((PC_Factor *)icc)->solvertype));
 82:   PetscCall(PCFactorClearComposedFunctions(pc));
 83:   PetscCall(PetscFree(pc->data));
 84:   PetscFunctionReturn(PETSC_SUCCESS);
 85: }

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

 91:   PetscFunctionBegin;
 92:   PetscCall(MatSolve(((PC_Factor *)icc)->fact, x, y));
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

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

100:   PetscFunctionBegin;
101:   PetscCall(MatMatSolve(((PC_Factor *)icc)->fact, X, Y));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

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

109:   PetscFunctionBegin;
110:   PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y));
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc, Vec x, Vec y)
115: {
116:   PC_ICC *icc = (PC_ICC *)pc->data;

118:   PetscFunctionBegin;
119:   PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y));
120:   PetscFunctionReturn(PETSC_SUCCESS);
121: }

123: static PetscErrorCode PCSetFromOptions_ICC(PC pc, PetscOptionItems *PetscOptionsObject)
124: {
125:   PC_ICC   *icc = (PC_ICC *)pc->data;
126:   PetscBool flg;
127:   /* PetscReal      dt[3];*/

129:   PetscFunctionBegin;
130:   PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options");
131:   PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));

133:   PetscCall(PetscOptionsReal("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", ((PC_Factor *)icc)->info.levels, &((PC_Factor *)icc)->info.levels, &flg));
134:   /*dt[0] = ((PC_Factor*)icc)->info.dt;
135:   dt[1] = ((PC_Factor*)icc)->info.dtcol;
136:   dt[2] = ((PC_Factor*)icc)->info.dtcount;
137:   PetscInt       dtmax = 3;
138:   PetscCall(PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg));
139:   if (flg) {
140:     PetscCall(PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]));
141:   }
142:   */
143:   PetscOptionsHeadEnd();
144:   PetscFunctionReturn(PETSC_SUCCESS);
145: }

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

149: /*MC
150:      PCICC - Incomplete Cholesky factorization preconditioners {cite}`chan1997approximate`

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

159:    Level: beginner

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

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

166:    By default, the Manteuffel shift {cite}`manteuffel1979shifted` is applied, for matrices with block size 1 only. Call `PCFactorSetShiftType`(pc,`MAT_SHIFT_POSITIVE_DEFINITE`);
167:    to turn off the shift.

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

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

179:   PetscFunctionBegin;
180:   PetscCall(PetscNew(&icc));
181:   pc->data = (void *)icc;
182:   PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
199: }