Actual source code: icc.c

petsc-master 2016-07-25
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;
 10:   PetscErrorCode         ierr;
 11:   MatInfo                info;
 12:   const MatSolverPackage stype;
 13:   MatFactorError         err;

 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:   MatFactorGetError(((PC_Factor*)icc)->fact,&err);
 36:   if (err) { /* FactorSymbolic() fails */
 37:     pc->failedreason = (PCFailedReason)err;
 38:     return(0);
 39:   }
 40: 
 41:   MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
 42:   MatFactorGetError(((PC_Factor*)icc)->fact,&err);
 43:   if (err) { /* FactorNumeric() fails */
 44:     pc->failedreason = (PCFailedReason)err;
 45:   }

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

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

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

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

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

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

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

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

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

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

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

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

129:   PetscOptionsHead(PetscOptionsObject,"ICC Options");
130:   PCSetFromOptions_Factor(PetscOptionsObject,pc);

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

148: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
149: {

153:   PCView_Factor(pc,viewer);
154:   return(0);
155: }

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

161: PetscErrorCode  PCFactorGetUseInPlace_ICC(PC pc,PetscBool *flg)
162: {
164:   *flg = PETSC_FALSE;
165:   return(0);
166: }

168: /*MC
169:      PCICC - Incomplete Cholesky factorization preconditioners.

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

178:    Level: beginner

180:   Concepts: incomplete Cholesky factorization

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

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

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

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

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


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

202: M*/

206: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
207: {
209:   PC_ICC         *icc;

212:   PetscNewLog(pc,&icc);

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

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

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

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

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

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