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:   PetscFunctionBegin;
 15:   pc->failedreason = PC_NOERROR;

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

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

 43:   PetscCall(ISDestroy(&cperm));
 44:   PetscCall(ISDestroy(&perm));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

131:   PetscFunctionBegin;
132:   PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options");
133:   PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));

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

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

151: /*MC
152:      PCICC - Incomplete Cholesky factorization preconditioners.

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

161:    Level: beginner

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

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

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

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

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

178: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCILU`, `PCLU`, `PCCHOLESKY`,
179:           `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
180:           `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
181:           `PCFactorSetLevels()`
182: M*/

184: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
185: {
186:   PC_ICC *icc;

188:   PetscFunctionBegin;
189:   PetscCall(PetscNew(&icc));
190:   pc->data = (void *)icc;
191:   PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ICC));

193:   ((PC_Factor *)icc)->info.fill      = 1.0;
194:   ((PC_Factor *)icc)->info.dtcol     = PETSC_DEFAULT;
195:   ((PC_Factor *)icc)->info.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;

197:   pc->ops->apply               = PCApply_ICC;
198:   pc->ops->matapply            = PCMatApply_ICC;
199:   pc->ops->applytranspose      = PCApply_ICC;
200:   pc->ops->setup               = PCSetUp_ICC;
201:   pc->ops->reset               = PCReset_ICC;
202:   pc->ops->destroy             = PCDestroy_ICC;
203:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
204:   pc->ops->view                = PCView_Factor;
205:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
206:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }