Actual source code: icc.c

petsc-main 2021-04-20
Report Typos and Errors

  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;
 12:   PetscBool              canuseordering;

 15:   pc->failedreason = PC_NOERROR;

 17:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
 18:   if (!pc->setupcalled) {
 19:     if (!((PC_Factor*)icc)->fact) {
 20:       MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
 21:     }
 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;

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

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

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

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

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

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

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

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

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

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

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

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

137:   PetscOptionsHead(PetscOptionsObject,"ICC Options");
138:   PCSetFromOptions_Factor(PetscOptionsObject,pc);

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

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

156: /*MC
157:      PCICC - Incomplete Cholesky factorization preconditioners.

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

166:    Level: beginner

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

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

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

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

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


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

189: M*/

191: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
192: {
194:   PC_ICC         *icc;

197:   PetscNewLog(pc,&icc);
198:   pc->data = (void*)icc;
199:   PCFactorInitialize(pc, 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: }