Actual source code: icc.c

petsc-master 2016-09-24
Report Typos and Errors
 2:  #include <../src/ksp/pc/impls/factor/icc/icc.h>

  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:   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->hdr.actualfill = info.fill_ratio_needed;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

160: /*MC
161:      PCICC - Incomplete Cholesky factorization preconditioners.

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

170:    Level: beginner

172:   Concepts: incomplete Cholesky factorization

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

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

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

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

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


189: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
190:            PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
191:            PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
192:            PCFactorSetLevels()

194: M*/

198: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
199: {
201:   PC_ICC         *icc;

204:   PetscNewLog(pc,&icc);
205:   pc->data = (void*)icc;
206:   PCFactorInitialize(pc);
207:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);

209:   ((PC_Factor*)icc)->factortype     = MAT_FACTOR_ICC;
210:   ((PC_Factor*)icc)->info.fill      = 1.0;
211:   ((PC_Factor*)icc)->info.dtcol     = PETSC_DEFAULT;
212:   ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;

214:   pc->ops->apply               = PCApply_ICC;
215:   pc->ops->applytranspose      = PCApply_ICC;
216:   pc->ops->setup               = PCSetUp_ICC;
217:   pc->ops->reset               = PCReset_ICC;
218:   pc->ops->destroy             = PCDestroy_ICC;
219:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
220:   pc->ops->view                = PCView_ICC;
221:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
222:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
223:   return(0);
224: }