Actual source code: icc.c

petsc-master 2017-07-23
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,cperm;
  8:   PetscErrorCode         ierr;
  9:   MatInfo                info;
 10:   const MatSolverPackage stype;
 11:   MatFactorError         err;

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

 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:     MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
 23:   } else if (pc->flag != SAME_NONZERO_PATTERN) {
 24:     MatDestroy(&((PC_Factor*)icc)->fact);
 25:     MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
 26:     MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
 27:   }
 28:   MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
 29:   icc->hdr.actualfill = info.fill_ratio_needed;

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

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

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

 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: }

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

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

 78: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
 79: {
 80:   PC_ICC         *icc = (PC_ICC*)pc->data;

 84:   MatSolve(((PC_Factor*)icc)->fact,x,y);
 85:   return(0);
 86: }

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

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

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

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

108: static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc)
109: {
110:   PC_ICC         *icc = (PC_ICC*)pc->data;
111:   PetscBool      flg;
113:   /* PetscReal      dt[3];*/

116:   PetscOptionsHead(PetscOptionsObject,"ICC Options");
117:   PCSetFromOptions_Factor(PetscOptionsObject,pc);

119:   PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
120:   /*dt[0] = ((PC_Factor*)icc)->info.dt;
121:   dt[1] = ((PC_Factor*)icc)->info.dtcol;
122:   dt[2] = ((PC_Factor*)icc)->info.dtcount;
123:   PetscInt       dtmax = 3;
124:   PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
125:   if (flg) {
126:     PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
127:   }
128:   */
129:   PetscOptionsTail();
130:   return(0);
131: }

133: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
134: {

138:   PCView_Factor(pc,viewer);
139:   return(0);
140: }

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

144: /*MC
145:      PCICC - Incomplete Cholesky factorization preconditioners.

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

154:    Level: beginner

156:   Concepts: incomplete Cholesky factorization

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

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

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

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

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


173: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
174:            PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
175:            PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
176:            PCFactorSetLevels()

178: M*/

180: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
181: {
183:   PC_ICC         *icc;

186:   PetscNewLog(pc,&icc);
187:   pc->data = (void*)icc;
188:   PCFactorInitialize(pc);
189:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);

191:   ((PC_Factor*)icc)->factortype     = MAT_FACTOR_ICC;
192:   ((PC_Factor*)icc)->info.fill      = 1.0;
193:   ((PC_Factor*)icc)->info.dtcol     = PETSC_DEFAULT;
194:   ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;

196:   pc->ops->apply               = PCApply_ICC;
197:   pc->ops->applytranspose      = PCApply_ICC;
198:   pc->ops->setup               = PCSetUp_ICC;
199:   pc->ops->reset               = PCReset_ICC;
200:   pc->ops->destroy             = PCDestroy_ICC;
201:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
202:   pc->ops->view                = PCView_ICC;
203:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
204:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
205:   return(0);
206: }