Actual source code: icc.c

  1: #define PETSCKSP_DLL

  3: /*
  4:    Defines a Cholesky factorization preconditioner for any Mat implementation.
  5:   Presently only provided for MPIRowbs format (i.e. BlockSolve).
  6: */

 8:  #include src/ksp/pc/impls/factor/icc/icc.h

 13: PetscErrorCode  PCFactorSetZeroPivot_ICC(PC pc,PetscReal z)
 14: {
 15:   PC_ICC *icc;

 18:   icc                 = (PC_ICC*)pc->data;
 19:   icc->info.zeropivot = z;
 20:   return(0);
 21: }

 27: PetscErrorCode  PCFactorSetShiftNonzero_ICC(PC pc,PetscReal shift)
 28: {
 29:   PC_ICC *dir;

 32:   dir = (PC_ICC*)pc->data;
 33:   if (shift == (PetscReal) PETSC_DECIDE) {
 34:     dir->info.shiftnz = 1.e-12;
 35:   } else {
 36:     dir->info.shiftnz = shift;
 37:   }
 38:   return(0);
 39: }

 45: PetscErrorCode  PCFactorSetShiftPd_ICC(PC pc,PetscTruth shift)
 46: {
 47:   PC_ICC *dir;
 48: 
 50:   dir = (PC_ICC*)pc->data;
 51:   if (shift) {
 52:     dir->info.shift_fraction = 0.0;
 53:     dir->info.shiftpd = 1.0;
 54:   } else {
 55:     dir->info.shiftpd = 0.0;
 56:   }
 57:   return(0);
 58: }

 64: PetscErrorCode  PCFactorSetMatOrdering_ICC(PC pc,MatOrderingType ordering)
 65: {
 66:   PC_ICC         *dir = (PC_ICC*)pc->data;
 68: 
 70:   PetscStrfree(dir->ordering);
 71:   PetscStrallocpy(ordering,&dir->ordering);
 72:   return(0);
 73: }

 79: PetscErrorCode  PCFactorSetFill_ICC(PC pc,PetscReal fill)
 80: {
 81:   PC_ICC *dir;

 84:   dir            = (PC_ICC*)pc->data;
 85:   dir->info.fill = fill;
 86:   return(0);
 87: }

 93: PetscErrorCode  PCFactorSetLevels_ICC(PC pc,PetscInt levels)
 94: {
 95:   PC_ICC *icc;

 98:   icc = (PC_ICC*)pc->data;
 99:   icc->info.levels = levels;
100:   return(0);
101: }

106: static PetscErrorCode PCSetup_ICC(PC pc)
107: {
108:   PC_ICC         *icc = (PC_ICC*)pc->data;
109:   IS             perm,cperm;
111:   MatInfo        info;

114:   MatGetOrdering(pc->pmat,icc->ordering,&perm,&cperm);

116:   if (!pc->setupcalled) {
117:     MatICCFactorSymbolic(pc->pmat,perm,&icc->info,&icc->fact);
118:   } else if (pc->flag != SAME_NONZERO_PATTERN) {
119:     MatDestroy(icc->fact);
120:     MatICCFactorSymbolic(pc->pmat,perm,&icc->info,&icc->fact);
121:   }
122:   MatGetInfo(icc->fact,MAT_LOCAL,&info);
123:   icc->actualfill = info.fill_ratio_needed;

125:   ISDestroy(cperm);
126:   ISDestroy(perm);
127:   MatCholeskyFactorNumeric(pc->pmat,&icc->info,&icc->fact);
128:   return(0);
129: }

133: static PetscErrorCode PCDestroy_ICC(PC pc)
134: {
135:   PC_ICC         *icc = (PC_ICC*)pc->data;

139:   if (icc->fact) {MatDestroy(icc->fact);}
140:   PetscStrfree(icc->ordering);
141:   PetscFree(icc);
142:   return(0);
143: }

147: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
148: {
149:   PC_ICC         *icc = (PC_ICC*)pc->data;

153:   MatSolve(icc->fact,x,y);
154:   return(0);
155: }

159: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
160: {
162:   PC_ICC         *icc = (PC_ICC*)pc->data;

165:   MatForwardSolve(icc->fact,x,y);
166:   return(0);
167: }

171: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
172: {
174:   PC_ICC         *icc = (PC_ICC*)pc->data;

177:   MatBackwardSolve(icc->fact,x,y);
178:   return(0);
179: }

183: static PetscErrorCode PCGetFactoredMatrix_ICC(PC pc,Mat *mat)
184: {
185:   PC_ICC *icc = (PC_ICC*)pc->data;

188:   *mat = icc->fact;
189:   return(0);
190: }

194: static PetscErrorCode PCSetFromOptions_ICC(PC pc)
195: {
196:   PC_ICC         *icc = (PC_ICC*)pc->data;
197:   char           tname[256];
198:   PetscTruth     flg;
200:   PetscFList     ordlist;

203:   MatOrderingRegisterAll(PETSC_NULL);
204:   PetscOptionsHead("ICC Options");
205:     PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",icc->info.levels,&icc->info.levels,&flg);
206:     PetscOptionsReal("-pc_factor_fill","Expected fill in factorization","PCFactorSetFill",icc->info.fill,&icc->info.fill,&flg);
207:     MatGetOrderingList(&ordlist);
208:     PetscOptionsList("-pc_factor_mat_ordering_type","Reorder to reduce nonzeros in ICC","PCFactorSetMatOrdering",ordlist,icc->ordering,tname,256,&flg);
209:     if (flg) {
210:       PCFactorSetMatOrdering(pc,tname);
211:     }
212:     PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
213:     if (flg) {
214:       PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);
215:     }
216:     PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",icc->info.shiftnz,&icc->info.shiftnz,0);
217:     PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShift",&flg);
218:     if (flg) {
219:       PCFactorSetShiftPd(pc,PETSC_TRUE);
220:     } else {
221:       PCFactorSetShiftPd(pc,PETSC_FALSE);
222:     }
223:     PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",icc->info.zeropivot,&icc->info.zeropivot,0);
224: 
225:   PetscOptionsTail();
226:   return(0);
227: }

231: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
232: {
233:   PC_ICC         *icc = (PC_ICC*)pc->data;
235:   PetscTruth     isstring,iascii;

238:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
239:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
240:   if (iascii) {
241:     if (icc->info.levels == 1) {
242:         PetscViewerASCIIPrintf(viewer,"  ICC: %D level of fill\n",(PetscInt)icc->info.levels);
243:     } else {
244:         PetscViewerASCIIPrintf(viewer,"  ICC: %D levels of fill\n",(PetscInt)icc->info.levels);
245:     }
246:     PetscViewerASCIIPrintf(viewer,"  ICC: factor fill ratio allocated %G\n",icc->info.fill);
247:     if (icc->info.shiftpd) {PetscViewerASCIIPrintf(viewer,"  ICC: using Manteuffel shift\n");}
248:     if (icc->fact) {
249:       PetscViewerASCIIPrintf(viewer,"  ICC: factor fill ratio needed %G\n",icc->actualfill);
250:       PetscViewerASCIIPrintf(viewer,"       Factored matrix follows\n");
251:       PetscViewerASCIIPushTab(viewer);
252:       PetscViewerASCIIPushTab(viewer);
253:       PetscViewerASCIIPushTab(viewer);
254:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
255:       MatView(icc->fact,viewer);
256:       PetscViewerPopFormat(viewer);
257:       PetscViewerASCIIPopTab(viewer);
258:       PetscViewerASCIIPopTab(viewer);
259:       PetscViewerASCIIPopTab(viewer);
260:     }
261:   } else if (isstring) {
262:     PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)icc->info.levels);
263:   } else {
264:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name);
265:   }
266:   return(0);
267: }

269: /*MC
270:      PCICC - Incomplete Cholesky factorization preconditioners.

272:    Options Database Keys:
273: +  -pc_factor_levels <k> - number of levels of fill for ICC(k)
274: .  -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
275:                       its factorization (overwrites original matrix)
276: .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
277: .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
278: .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
279: -  -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
280:    is optional with PETSC_TRUE being the default

282:    Level: beginner

284:   Concepts: incomplete Cholesky factorization

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

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

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

292:           By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftPd(pc,PETSC_FALSE);
293:           to turn off the shift.


296: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
297:            PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(), 
298:            PCFactorSetFill(), PCFactorSetMatOrdering(), PCFactorSetReuseOrdering(), 
299:            PCFactorSetLevels(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd(),

301: M*/

306: PetscErrorCode  PCCreate_ICC(PC pc)
307: {
309:   PC_ICC         *icc;

312:   PetscNew(PC_ICC,&icc);
313:   PetscLogObjectMemory(pc,sizeof(PC_ICC));

315:   icc->fact                  = 0;
316:   PetscStrallocpy(MATORDERING_NATURAL,&icc->ordering);
317:   MatFactorInfoInitialize(&icc->info);
318:   icc->info.levels          = 0;
319:   icc->info.fill          = 1.0;
320:   icc->implctx            = 0;

322:   icc->info.dtcol              = PETSC_DEFAULT;
323:   icc->info.shiftnz            = 0.0;
324:   icc->info.shiftpd            = 1.0; /* true */
325:   icc->info.shift_fraction     = 0.0;
326:   icc->info.zeropivot          = 1.e-12;
327:   pc->data                       = (void*)icc;

329:   pc->ops->apply               = PCApply_ICC;
330:   pc->ops->setup               = PCSetup_ICC;
331:   pc->ops->destroy               = PCDestroy_ICC;
332:   pc->ops->setfromoptions      = PCSetFromOptions_ICC;
333:   pc->ops->view                = PCView_ICC;
334:   pc->ops->getfactoredmatrix   = PCGetFactoredMatrix_ICC;
335:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
336:   pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;

338:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_ICC",
339:                     PCFactorSetZeroPivot_ICC);
340:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_ICC",
341:                     PCFactorSetShiftNonzero_ICC);
342:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_ICC",
343:                     PCFactorSetShiftPd_ICC);

345:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_ICC",
346:                     PCFactorSetLevels_ICC);
347:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_ICC",
348:                     PCFactorSetFill_ICC);
349:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrdering_C","PCFactorSetMatOrdering_ICC",
350:                     PCFactorSetMatOrdering_ICC);
351:   return(0);
352: }