Actual source code: cholesky.c

  1: #define PETSCKSP_DLL

  3: /*
  4:    Defines a direct factorization preconditioner for any Mat implementation
  5:    Note: this need not be consided a preconditioner since it supplies
  6:          a direct solver.
  7: */
 8:  #include private/pcimpl.h

 10: typedef struct {
 11:   Mat             fact;             /* factored matrix */
 12:   PetscReal       actualfill;       /* actual fill in factor */
 13:   PetscTruth      inplace;          /* flag indicating in-place factorization */
 14:   IS              row,col;          /* index sets used for reordering */
 15:   MatOrderingType ordering;         /* matrix ordering */
 16:   PetscTruth      reuseordering;    /* reuses previous reordering computed */
 17:   PetscTruth      reusefill;        /* reuse fill from previous Cholesky */
 18:   MatFactorInfo   info;
 19: } PC_Cholesky;

 24: PetscErrorCode  PCFactorSetZeroPivot_Cholesky(PC pc,PetscReal z)
 25: {
 26:   PC_Cholesky *ch;

 29:   ch                 = (PC_Cholesky*)pc->data;
 30:   ch->info.zeropivot = z;
 31:   return(0);
 32: }

 38: PetscErrorCode  PCFactorSetShiftNonzero_Cholesky(PC pc,PetscReal shift)
 39: {
 40:   PC_Cholesky *dir;

 43:   dir = (PC_Cholesky*)pc->data;
 44:   if (shift == (PetscReal) PETSC_DECIDE) {
 45:     dir->info.shiftnz = 1.e-12;
 46:   } else {
 47:     dir->info.shiftnz = shift;
 48:   }
 49:   return(0);
 50: }

 56: PetscErrorCode  PCFactorSetShiftPd_Cholesky(PC pc,PetscTruth shift)
 57: {
 58:   PC_Cholesky *dir;
 59: 
 61:   dir = (PC_Cholesky*)pc->data;
 62:   if (shift) {
 63:     dir->info.shift_fraction = 0.0;
 64:     dir->info.shiftpd = 1.0;
 65:   } else {
 66:     dir->info.shiftpd = 0.0;
 67:   }
 68:   return(0);
 69: }

 75: PetscErrorCode  PCFactorSetReuseOrdering_Cholesky(PC pc,PetscTruth flag)
 76: {
 77:   PC_Cholesky *lu;
 78: 
 80:   lu               = (PC_Cholesky*)pc->data;
 81:   lu->reuseordering = flag;
 82:   return(0);
 83: }

 89: PetscErrorCode  PCFactorSetReuseFill_Cholesky(PC pc,PetscTruth flag)
 90: {
 91:   PC_Cholesky *lu;
 92: 
 94:   lu = (PC_Cholesky*)pc->data;
 95:   lu->reusefill = flag;
 96:   return(0);
 97: }

102: static PetscErrorCode PCSetFromOptions_Cholesky(PC pc)
103: {
104:   PC_Cholesky    *lu = (PC_Cholesky*)pc->data;
106:   PetscTruth     flg;
107:   char           tname[256];
108:   PetscFList     ordlist;
109: 
111:   MatOrderingRegisterAll(PETSC_NULL);
112:   PetscOptionsHead("Cholesky options");
113:   PetscOptionsName("-pc_factor_in_place","Form Cholesky in the same memory as the matrix","PCFactorSetUseInPlace",&flg);
114:   if (flg) {
115:     PCFactorSetUseInPlace(pc);
116:   }
117:   PetscOptionsReal("-pc_factor_fill","Expected non-zeros in Cholesky/non-zeros in matrix","PCFactorSetFill",lu->info.fill,&lu->info.fill,0);
118: 
119:   PetscOptionsName("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",&flg);
120:   if (flg) {
121:     PCFactorSetReuseFill(pc,PETSC_TRUE);
122:   }
123:   PetscOptionsName("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",&flg);
124:   if (flg) {
125:     PCFactorSetReuseOrdering(pc,PETSC_TRUE);
126:   }
127: 
128:   MatGetOrderingList(&ordlist);
129:   PetscOptionsList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in Cholesky","PCFactorSetMatOrdering",ordlist,lu->ordering,tname,256,&flg);
130:   if (flg) {
131:     PCFactorSetMatOrdering(pc,tname);
132:   }
133:   PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
134:   if (flg) {
135:     PCFactorSetShiftNonzero(pc,(PetscReal) PETSC_DECIDE);
136:   }
137:   PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",lu->info.shiftnz,&lu->info.shiftnz,0);
138:   PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",&flg);
139:   if (flg) {
140:     PCFactorSetShiftPd(pc,PETSC_TRUE);
141:   }
142:   PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",lu->info.zeropivot,&lu->info.zeropivot,0);

144:   PetscOptionsTail();
145:   return(0);
146: }

150: static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer)
151: {
152:   PC_Cholesky    *lu = (PC_Cholesky*)pc->data;
154:   PetscTruth     iascii,isstring;
155: 
157:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
158:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
159:   if (iascii) {
160: 
161:     if (lu->inplace) {PetscViewerASCIIPrintf(viewer,"  Cholesky: in-place factorization\n");}
162:     else             {PetscViewerASCIIPrintf(viewer,"  Cholesky: out-of-place factorization\n");}
163:     PetscViewerASCIIPrintf(viewer,"    matrix ordering: %s\n",lu->ordering);
164:     if (lu->fact) {
165:       PetscViewerASCIIPrintf(viewer,"  Cholesky: factor fill ratio needed %G\n",lu->actualfill);
166:       PetscViewerASCIIPrintf(viewer,"       Factored matrix follows\n");
167:       PetscViewerASCIIPushTab(viewer);
168:       PetscViewerASCIIPushTab(viewer);
169:       PetscViewerASCIIPushTab(viewer);
170:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
171:       MatView(lu->fact,viewer);
172:       PetscViewerPopFormat(viewer);
173:       PetscViewerASCIIPopTab(viewer);
174:       PetscViewerASCIIPopTab(viewer);
175:       PetscViewerASCIIPopTab(viewer);
176:     }
177:     if (lu->reusefill)    {PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorization\n");}
178:     if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorization\n");}
179:   } else if (isstring) {
180:     PetscViewerStringSPrintf(viewer," order=%s",lu->ordering);
181:   } else {
182:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCCHOLESKY",((PetscObject)viewer)->type_name);
183:   }
184:   return(0);
185: }

189: static PetscErrorCode PCGetFactoredMatrix_Cholesky(PC pc,Mat *mat)
190: {
191:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;
192: 
194:   if (!dir->fact) SETERRQ(PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
195:   *mat = dir->fact;
196:   return(0);
197: }

201: static PetscErrorCode PCSetUp_Cholesky(PC pc)
202: {
204:   PetscTruth     flg;
205:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

208:   if (dir->reusefill && pc->setupcalled) dir->info.fill = dir->actualfill;
209: 
210:   if (dir->inplace) {
211:     if (dir->row && dir->col && (dir->row != dir->col)) {
212:       ISDestroy(dir->row);
213:       dir->row = 0;
214:     }
215:     if (dir->col) {
216:       ISDestroy(dir->col);
217:       dir->col = 0;
218:     }
219:     MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
220:     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
221:       ISDestroy(dir->col);
222:       dir->col=0;
223:     }
224:     if (dir->row) {PetscLogObjectParent(pc,dir->row);}
225:     MatCholeskyFactor(pc->pmat,dir->row,&dir->info);
226:     dir->fact = pc->pmat;
227:   } else {
228:     MatInfo info;
229:     if (!pc->setupcalled) {
230:       MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
231:       /* check if dir->row == dir->col */
232:       ISEqual(dir->row,dir->col,&flg);
233:       if (!flg) SETERRQ(PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
234:       ISDestroy(dir->col); /* only pass one ordering into CholeskyFactor */
235:       dir->col=0;

237:       PetscOptionsHasName(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);
238:       if (flg) {
239:         PetscReal tol = 1.e-10;
240:         PetscOptionsGetReal(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);
241:         MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
242:       }
243:       if (dir->row) {PetscLogObjectParent(pc,dir->row);}
244:       MatCholeskyFactorSymbolic(pc->pmat,dir->row,&dir->info,&dir->fact);
245:       MatGetInfo(dir->fact,MAT_LOCAL,&info);
246:       dir->actualfill = info.fill_ratio_needed;
247:       PetscLogObjectParent(pc,dir->fact);
248:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
249:       if (!dir->reuseordering) {
250:         if (dir->row && dir->col && (dir->row != dir->col)) {
251:           ISDestroy(dir->row);
252:           dir->row = 0;
253:         }
254:         if (dir->col) {
255:           ISDestroy(dir->col);
256:           dir->col =0;
257:         }
258:         MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
259:         if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
260:           ISDestroy(dir->col);
261:           dir->col=0;
262:         }
263:         PetscOptionsHasName(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);
264:         if (flg) {
265:           PetscReal tol = 1.e-10;
266:           PetscOptionsGetReal(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);
267:           MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
268:         }
269:         if (dir->row) {PetscLogObjectParent(pc,dir->row);}
270:       }
271:       MatDestroy(dir->fact);
272:       MatCholeskyFactorSymbolic(pc->pmat,dir->row,&dir->info,&dir->fact);
273:       MatGetInfo(dir->fact,MAT_LOCAL,&info);
274:       dir->actualfill = info.fill_ratio_needed;
275:       PetscLogObjectParent(pc,dir->fact);
276:     }
277:     MatCholeskyFactorNumeric(pc->pmat,&dir->info,&dir->fact);
278:   }
279:   return(0);
280: }

284: static PetscErrorCode PCDestroy_Cholesky(PC pc)
285: {
286:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

290:   if (!dir->inplace && dir->fact) {MatDestroy(dir->fact);}
291:   if (dir->row) {ISDestroy(dir->row);}
292:   if (dir->col) {ISDestroy(dir->col);}
293:   PetscStrfree(dir->ordering);
294:   PetscFree(dir);
295:   return(0);
296: }

300: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
301: {
302:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
304: 
306:   if (dir->inplace) {MatSolve(pc->pmat,x,y);}
307:   else              {MatSolve(dir->fact,x,y);}
308:   return(0);
309: }

313: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
314: {
315:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

319:   if (dir->inplace) {MatSolveTranspose(pc->pmat,x,y);}
320:   else              {MatSolveTranspose(dir->fact,x,y);}
321:   return(0);
322: }

324: /* -----------------------------------------------------------------------------------*/

329: PetscErrorCode  PCFactorSetFill_Cholesky(PC pc,PetscReal fill)
330: {
331:   PC_Cholesky *dir;
332: 
334:   dir = (PC_Cholesky*)pc->data;
335:   dir->info.fill = fill;
336:   return(0);
337: }

343: PetscErrorCode  PCFactorSetUseInPlace_Cholesky(PC pc)
344: {
345:   PC_Cholesky *dir;

348:   dir = (PC_Cholesky*)pc->data;
349:   dir->inplace = PETSC_TRUE;
350:   return(0);
351: }

357: PetscErrorCode  PCFactorSetMatOrdering_Cholesky(PC pc,MatOrderingType ordering)
358: {
359:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

363:   PetscStrfree(dir->ordering);
364:   PetscStrallocpy(ordering,&dir->ordering);
365:   return(0);
366: }

369: /* -----------------------------------------------------------------------------------*/

373: /*@
374:    PCFactorSetReuseOrdering - When similar matrices are factored, this
375:    causes the ordering computed in the first factor to be used for all
376:    following factors.

378:    Collective on PC

380:    Input Parameters:
381: +  pc - the preconditioner context
382: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

384:    Options Database Key:
385: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

387:    Level: intermediate

389: .keywords: PC, levels, reordering, factorization, incomplete, LU

391: .seealso: PCFactorSetReuseFill()
392: @*/
393: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscTruth flag)
394: {
395:   PetscErrorCode ierr,(*f)(PC,PetscTruth);

399:   PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",(void (**)(void))&f);
400:   if (f) {
401:     (*f)(pc,flag);
402:   }
403:   return(0);
404: }

408: /*@
409:    PCFactorSetReuseFill - When matrices with same different nonzero structure are factored,
410:    this causes later ones to use the fill ratio computed in the initial factorization.

412:    Collective on PC

414:    Input Parameters:
415: +  pc - the preconditioner context
416: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

418:    Options Database Key:
419: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()

421:    Level: intermediate

423: .keywords: PC, levels, reordering, factorization, incomplete, Cholesky

425: .seealso: PCFactorSetReuseOrdering()
426: @*/
427: PetscErrorCode  PCFactorSetReuseFill(PC pc,PetscTruth flag)
428: {
429:   PetscErrorCode ierr,(*f)(PC,PetscTruth);

433:   PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseFill_C",(void (**)(void))&f);
434:   if (f) {
435:     (*f)(pc,flag);
436:   }
437:   return(0);
438: }

440: /*MC
441:    PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner

443:    Options Database Keys:
444: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
445: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
446: .  -pc_factor_fill <fill> - Sets fill amount
447: .  -pc_factor_in_place - Activates in-place factorization
448: .  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
449: .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
450: -  -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
451:    is optional with PETSC_TRUE being the default

453:    Notes: Not all options work for all matrix formats

455:    Level: beginner

457:    Concepts: Cholesky factorization, direct solver

459:    Notes: Usually this will compute an "exact" solution in one iteration and does 
460:           not need a Krylov method (i.e. you can use -ksp_type preonly, or 
461:           KSPSetType(ksp,KSPPREONLY) for the Krylov method

463: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
464:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCGetFactoredMatrix(),
465:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(),
466:            PCFactorSetUseInPlace(), PCFactorSetMatOrdering(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd()

468: M*/

473: PetscErrorCode  PCCreate_Cholesky(PC pc)
474: {
476:   PC_Cholesky    *dir;

479:   PetscNew(PC_Cholesky,&dir);
480:   PetscLogObjectMemory(pc,sizeof(PC_Cholesky));

482:   dir->fact                   = 0;
483:   dir->inplace                = PETSC_FALSE;
484:   MatFactorInfoInitialize(&dir->info);
485:   dir->info.fill              = 5.0;
486:   dir->info.shiftnz           = 0.0;
487:   dir->info.shiftpd           = 0.0; /* false */
488:   dir->info.shift_fraction    = 0.0;
489:   dir->info.pivotinblocks     = 1.0;
490:   dir->col                    = 0;
491:   dir->row                    = 0;
492:   PetscStrallocpy(MATORDERING_NATURAL,&dir->ordering);
493:   dir->reusefill        = PETSC_FALSE;
494:   dir->reuseordering    = PETSC_FALSE;
495:   pc->data              = (void*)dir;

497:   pc->ops->destroy           = PCDestroy_Cholesky;
498:   pc->ops->apply             = PCApply_Cholesky;
499:   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
500:   pc->ops->setup             = PCSetUp_Cholesky;
501:   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
502:   pc->ops->view              = PCView_Cholesky;
503:   pc->ops->applyrichardson   = 0;
504:   pc->ops->getfactoredmatrix = PCGetFactoredMatrix_Cholesky;

506:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Cholesky",
507:                     PCFactorSetZeroPivot_Cholesky);
508:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Cholesky",
509:                     PCFactorSetShiftNonzero_Cholesky);
510:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Cholesky",
511:                     PCFactorSetShiftPd_Cholesky);

513:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Cholesky",
514:                     PCFactorSetFill_Cholesky);
515:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky",
516:                     PCFactorSetUseInPlace_Cholesky);
517:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrdering_C","PCFactorSetMatOrdering_Cholesky",
518:                     PCFactorSetMatOrdering_Cholesky);
519:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky",
520:                     PCFactorSetReuseOrdering_Cholesky);
521:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky",
522:                     PCFactorSetReuseFill_Cholesky);
523:   return(0);
524: }