Actual source code: icc.c

petsc-3.7.6 2017-04-24
Report Typos and Errors
  2: #include <../src/ksp/pc/impls/factor/icc/icc.h>   /*I "petscpc.h" I*/

  6: static PetscErrorCode PCSetUp_ICC(PC pc)
  7: {
  8:   PC_ICC         *icc = (PC_ICC*)pc->data;
  9:   IS             perm,cperm;
 11:   MatInfo        info;
 12:   Mat            F;
 13:   const MatSolverPackage stype;

 16:   pc->failedreason = PC_NOERROR;
 17:   MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);

 19:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
 20:   if (!pc->setupcalled) {
 21:     if (!((PC_Factor*)icc)->fact) {
 22:       MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
 23:     }
 24:     MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
 25:   } else if (pc->flag != SAME_NONZERO_PATTERN) {
 26:     MatDestroy(&((PC_Factor*)icc)->fact);
 27:     MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
 28:     MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
 29:   }
 30:   MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
 31:   icc->actualfill = info.fill_ratio_needed;

 33:   ISDestroy(&cperm);
 34:   ISDestroy(&perm);

 36:   F = ((PC_Factor*)icc)->fact;
 37:   if (F->errortype) { /* FactorSymbolic() fails */
 38:     pc->failedreason = (PCFailedReason)F->errortype;
 39:     return(0);
 40:   }
 41: 
 42:   MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
 43:   if (F->errortype) { /* FactorNumeric() fails */
 44:     pc->failedreason = (PCFailedReason)F->errortype;
 45:   }

 47:   PCFactorGetMatSolverPackage(pc,&stype);
 48:   if (!stype) {
 49:     PCFactorSetMatSolverPackage(pc,((PC_Factor*)icc)->fact->solvertype);
 50:   }
 51:   return(0);
 52: }

 56: static PetscErrorCode PCReset_ICC(PC pc)
 57: {
 58:   PC_ICC         *icc = (PC_ICC*)pc->data;

 62:   MatDestroy(&((PC_Factor*)icc)->fact);
 63:   return(0);
 64: }

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

 74:   PCReset_ICC(pc);
 75:   PetscFree(((PC_Factor*)icc)->ordering);
 76:   PetscFree(((PC_Factor*)icc)->solvertype);
 77:   PetscFree(pc->data);
 78:   return(0);
 79: }

 83: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
 84: {
 85:   PC_ICC         *icc = (PC_ICC*)pc->data;

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

 95: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
 96: {
 98:   PC_ICC         *icc = (PC_ICC*)pc->data;

101:   MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
102:   return(0);
103: }

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

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

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

127:   PetscOptionsHead(PetscOptionsObject,"ICC Options");
128:   PCSetFromOptions_Factor(PetscOptionsObject,pc);

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

146: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
147: {

151:   PCView_Factor(pc,viewer);
152:   return(0);
153: }

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

159: PetscErrorCode  PCFactorGetUseInPlace_ICC(PC pc,PetscBool *flg)
160: {
162:   *flg = PETSC_FALSE;
163:   return(0);
164: }

166: /*MC
167:      PCICC - Incomplete Cholesky factorization preconditioners.

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

176:    Level: beginner

178:   Concepts: incomplete Cholesky factorization

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

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

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

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

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


195: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
196:            PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
197:            PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
198:            PCFactorSetLevels()

200: M*/

204: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
205: {
207:   PC_ICC         *icc;

210:   PetscNewLog(pc,&icc);

212:   ((PC_Factor*)icc)->fact = 0;

214:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);
215:   MatFactorInfoInitialize(&((PC_Factor*)icc)->info);

217:   ((PC_Factor*)icc)->factortype  = MAT_FACTOR_ICC;
218:   ((PC_Factor*)icc)->info.levels = 0.;
219:   ((PC_Factor*)icc)->info.fill   = 1.0;
220:   icc->implctx                   = 0;

222:   ((PC_Factor*)icc)->info.dtcol       = PETSC_DEFAULT;
223:   ((PC_Factor*)icc)->info.shifttype   = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
224:   ((PC_Factor*)icc)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
225:   ((PC_Factor*)icc)->info.zeropivot   = 100.0*PETSC_MACHINE_EPSILON;

227:   pc->data                     = (void*)icc;
228:   pc->ops->apply               = PCApply_ICC;
229:   pc->ops->applytranspose      = PCApply_ICC;
230:   pc->ops->setup               = PCSetUp_ICC;
231:   pc->ops->reset               = PCReset_ICC;
232:   pc->ops->destroy             = PCDestroy_ICC;
233:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
234:   pc->ops->view                = PCView_ICC;
235:   pc->ops->getfactoredmatrix   = PCFactorGetMatrix_Factor;
236:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
237:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;

239:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
240:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetZeroPivot_C",PCFactorGetZeroPivot_Factor);
241:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
242:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftType_C",PCFactorGetShiftType_Factor);
243:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
244:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftAmount_C",PCFactorGetShiftAmount_Factor);
245:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
246:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
247:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C",PCFactorSetLevels_Factor);
248:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetLevels_C",PCFactorGetLevels_Factor);
249:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
250:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
251:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
252:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);
253:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_ICC);
254:   return(0);
255: }