Actual source code: icc.c

petsc-3.3-p7 2013-05-11
  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;

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

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

 29:   ISDestroy(&cperm);
 30:   ISDestroy(&perm);
 31:   MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
 32:   return(0);
 33: }

 37: static PetscErrorCode PCReset_ICC(PC pc)
 38: {
 39:   PC_ICC         *icc = (PC_ICC*)pc->data;

 43:   MatDestroy(&((PC_Factor*)icc)->fact);
 44:   return(0);
 45: }

 49: static PetscErrorCode PCDestroy_ICC(PC pc)
 50: {
 51:   PC_ICC         *icc = (PC_ICC*)pc->data;

 55:   PCReset_ICC(pc);
 56:   PetscFree(((PC_Factor*)icc)->ordering);
 57:   PetscFree(((PC_Factor*)icc)->solvertype);
 58:   PetscFree(pc->data);
 59:   return(0);
 60: }

 64: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
 65: {
 66:   PC_ICC         *icc = (PC_ICC*)pc->data;

 70:   MatSolve(((PC_Factor*)icc)->fact,x,y);
 71:   return(0);
 72: }

 76: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
 77: {
 79:   PC_ICC         *icc = (PC_ICC*)pc->data;

 82:   MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
 83:   return(0);
 84: }

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

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

100: static PetscErrorCode PCSetFromOptions_ICC(PC pc)
101: {
102:   PC_ICC         *icc = (PC_ICC*)pc->data;
103:   PetscBool      flg;
104:   /* PetscReal      dt[3];*/

108:   PetscOptionsHead("ICC Options");
109:     PCSetFromOptions_Factor(pc);

111:     PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
112:     /*dt[0] = ((PC_Factor*)icc)->info.dt;
113:     dt[1] = ((PC_Factor*)icc)->info.dtcol;
114:     dt[2] = ((PC_Factor*)icc)->info.dtcount;
115:     PetscInt       dtmax = 3;
116:     PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
117:     if (flg) {
118:       PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
119:     }
120:     */
121:   PetscOptionsTail();
122:   return(0);
123: }

127: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
128: {
130: 
132:   PCView_Factor(pc,viewer);
133:   return(0);
134: }

136: EXTERN_C_BEGIN
137: extern PetscErrorCode  PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);
138: EXTERN_C_END

140: /*MC
141:      PCICC - Incomplete Cholesky factorization preconditioners.

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

150:    Level: beginner

152:   Concepts: incomplete Cholesky factorization

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

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

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

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

163:    References:
164:    Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST
165:       http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical
166:       Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
167:       Science and Engineering, Kluwer, pp. 167--202.


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

175: M*/

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

186:   PetscNewLog(pc,PC_ICC,&icc);

188:   ((PC_Factor*)icc)->fact                  = 0;
189:   PetscStrallocpy(MATORDERINGNATURAL,&((PC_Factor*)icc)->ordering);
190:   PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)icc)->solvertype);
191:   MatFactorInfoInitialize(&((PC_Factor*)icc)->info);
192:   ((PC_Factor*)icc)->factortype         = MAT_FACTOR_ICC;
193:   ((PC_Factor*)icc)->info.levels        = 0.;
194:   ((PC_Factor*)icc)->info.fill          = 1.0;
195:   icc->implctx            = 0;

197:   ((PC_Factor*)icc)->info.dtcol       = PETSC_DEFAULT;
198:   ((PC_Factor*)icc)->info.shifttype   = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
199:   ((PC_Factor*)icc)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
200:   ((PC_Factor*)icc)->info.zeropivot   = 100.0*PETSC_MACHINE_EPSILON;

202:   pc->data                       = (void*)icc;
203:   pc->ops->apply               = PCApply_ICC;
204:   pc->ops->applytranspose      = PCApply_ICC;
205:   pc->ops->setup               = PCSetup_ICC;
206:   pc->ops->reset                 = PCReset_ICC;
207:   pc->ops->destroy               = PCDestroy_ICC;
208:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
209:   pc->ops->view                = PCView_ICC;
210:   pc->ops->getfactoredmatrix   = PCFactorGetMatrix_Factor;
211:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
212:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;

214:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor",
215:                     PCFactorSetUpMatSolverPackage_Factor);
216:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",
217:                     PCFactorGetMatSolverPackage_Factor);
218:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",
219:                     PCFactorSetZeroPivot_Factor);
220:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor",
221:                     PCFactorSetShiftType_Factor);
222:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor",
223:                     PCFactorSetShiftAmount_Factor);
224:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor",
225:                     PCFactorSetLevels_Factor);
226:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",
227:                     PCFactorSetFill_Factor);
228:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",
229:                     PCFactorSetMatOrderingType_Factor);
230:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",
231:                     PCFactorSetMatSolverPackage_Factor);
232:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetDropTolerance_C","PCFactorSetDropTolerance_ILU",
233:                     PCFactorSetDropTolerance_ILU);
234:   return(0);
235: }
236: EXTERN_C_END