Actual source code: icc.c

petsc-3.7.3 2016-07-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:   MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

175:    Level: beginner

177:   Concepts: incomplete Cholesky factorization

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

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

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

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

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


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

199: M*/

203: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
204: {
206:   PC_ICC         *icc;

209:   PetscNewLog(pc,&icc);

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

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

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

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

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

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