Actual source code: cholesky.c

  1: /*$Id: cholesky.c,v 1.12 2001/08/06 21:16:36 bsmith Exp $*/
  2: /*
  3:    Defines a direct factorization preconditioner for any Mat implementation
  4:    Note: this need not be consided a preconditioner since it supplies
  5:          a direct solver.
  6: */
 7:  #include src/sles/pc/pcimpl.h

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

 20: EXTERN_C_BEGIN
 21: #undef __FUNCT__  
 23: int PCCholeskySetReuseOrdering_Cholesky(PC pc,PetscTruth flag)
 24: {
 25:   PC_Cholesky *lu;
 26: 
 28:   lu               = (PC_Cholesky*)pc->data;
 29:   lu->reuseordering = flag;
 30:   return(0);
 31: }
 32: EXTERN_C_END

 34: EXTERN_C_BEGIN
 35: #undef __FUNCT__  
 37: int PCCholeskySetReuseFill_Cholesky(PC pc,PetscTruth flag)
 38: {
 39:   PC_Cholesky *lu;
 40: 
 42:   lu = (PC_Cholesky*)pc->data;
 43:   lu->reusefill = flag;
 44:   return(0);
 45: }
 46: EXTERN_C_END

 48: #undef __FUNCT__  
 50: static int PCSetFromOptions_Cholesky(PC pc)
 51: {
 52:   PC_Cholesky *lu = (PC_Cholesky*)pc->data;
 53:   int         ierr;
 54:   PetscTruth  flg;
 55:   char        tname[256];
 56:   PetscFList  ordlist;
 57: 
 59:   MatOrderingRegisterAll(PETSC_NULL);
 60:   PetscOptionsHead("Cholesky options");
 61:   PetscOptionsName("-pc_cholesky_in_place","Form Cholesky in the same memory as the matrix","PCCholeskySetUseInPlace",&flg);
 62:   if (flg) {
 63:     PCCholeskySetUseInPlace(pc);
 64:   }
 65:   PetscOptionsReal("-pc_cholesky_fill","Expected non-zeros in Cholesky/non-zeros in matrix","PCCholeskySetFill",lu->info.fill,&lu->info.fill,0);
 66: 
 67:   PetscOptionsName("-pc_cholesky_reuse_fill","Use fill from previous factorization","PCCholeskySetReuseFill",&flg);
 68:   if (flg) {
 69:     PCCholeskySetReuseFill(pc,PETSC_TRUE);
 70:   }
 71:   PetscOptionsName("-pc_cholesky_reuse_ordering","Reuse ordering from previous factorization","PCCholeskySetReuseOrdering",&flg);
 72:   if (flg) {
 73:     PCCholeskySetReuseOrdering(pc,PETSC_TRUE);
 74:   }
 75: 
 76:   MatGetOrderingList(&ordlist);
 77:   PetscOptionsList("-pc_cholesky_mat_ordering_type","Reordering to reduce nonzeros in Cholesky","PCCholeskySetMatOrdering",ordlist,lu->ordering,tname,256,&flg);
 78:   if (flg) {
 79:     PCCholeskySetMatOrdering(pc,tname);
 80:   }
 81: 
 82:   PetscOptionsTail();
 83:   return(0);
 84: }

 86: #undef __FUNCT__  
 88: static int PCView_Cholesky(PC pc,PetscViewer viewer)
 89: {
 90:   PC_Cholesky *lu = (PC_Cholesky*)pc->data;
 91:   int         ierr;
 92:   PetscTruth  isascii,isstring;
 93: 
 95:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
 96:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
 97:   if (isascii) {
 98:     MatInfo info;
 99: 
100:     if (lu->inplace) {PetscViewerASCIIPrintf(viewer,"  Cholesky: in-place factorizationn");}
101:     else             {PetscViewerASCIIPrintf(viewer,"  Cholesky: out-of-place factorizationn");}
102:     PetscViewerASCIIPrintf(viewer,"    matrix ordering: %sn",lu->ordering);
103:     if (lu->fact) {
104:       MatGetInfo(lu->fact,MAT_LOCAL,&info);
105:       PetscViewerASCIIPrintf(viewer,"    Cholesky nonzeros %gn",info.nz_used);
106:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_FACTOR_INFO);
107:       MatView(lu->fact,viewer);
108:       PetscViewerPopFormat(viewer);
109:     }
110:     if (lu->reusefill)    {PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorizationn");}
111:     if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorizationn");}
112:   } else if (isstring) {
113:     PetscViewerStringSPrintf(viewer," order=%s",lu->ordering);
114:   } else {
115:     SETERRQ1(1,"Viewer type %s not supported for PCCholesky",((PetscObject)viewer)->type_name);
116:   }
117:   return(0);
118: }

120: #undef __FUNCT__  
122: static int PCGetFactoredMatrix_Cholesky(PC pc,Mat *mat)
123: {
124:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;
125: 
127:   if (!dir->fact) SETERRQ(1,"Matrix not yet factored; call after SLESSetUp() or PCSetUp()");
128:   *mat = dir->fact;
129:   return(0);
130: }

132: #undef __FUNCT__  
134: static int PCSetUp_Cholesky(PC pc)
135: {
136:   int        ierr;
137:   PetscTruth flg;
138:   PC_Cholesky      *dir = (PC_Cholesky*)pc->data;

141:   if (dir->reusefill && pc->setupcalled) dir->info.fill = dir->actualfill;
142: 
143:   if (dir->inplace) {
144:     if (dir->row && dir->col && (dir->row != dir->col)) {
145:       ISDestroy(dir->row);
146:       dir->row = 0;
147:     }
148:     if (dir->col) {
149:       ISDestroy(dir->col);
150:       dir->col = 0;
151:     }
152:     MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
153:     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
154:       ISDestroy(dir->col);
155:       dir->col=0;
156:     }
157:     if (dir->row) {PetscLogObjectParent(pc,dir->row);}
158:     MatCholeskyFactor(pc->pmat,dir->row,dir->info.fill);
159:     dir->fact = pc->pmat;
160:   } else {
161:     MatInfo info;
162:     if (!pc->setupcalled) {
163:       MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
164:       if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
165:         ISDestroy(dir->col);
166:         dir->col=0;
167:       }
168:       PetscOptionsHasName(pc->prefix,"-pc_cholesky_nonzeros_along_diagonal",&flg);
169:       if (flg) {
170:         PetscReal tol = 1.e-10;
171:         PetscOptionsGetReal(pc->prefix,"-pc_cholesky_nonzeros_along_diagonal",&tol,PETSC_NULL);
172:         MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
173:       }
174:       if (dir->row) {PetscLogObjectParent(pc,dir->row);}
175:       MatCholeskyFactorSymbolic(pc->pmat,dir->row,dir->info.fill,&dir->fact);
176:       MatGetInfo(dir->fact,MAT_LOCAL,&info);
177:       dir->actualfill = info.fill_ratio_needed;
178:       PetscLogObjectParent(pc,dir->fact);
179:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
180:       if (!dir->reuseordering) {
181:         if (dir->row && dir->col && (dir->row != dir->col)) {
182:           ISDestroy(dir->row);
183:           dir->row = 0;
184:         }
185:         if (dir->col) {
186:           ISDestroy(dir->col);
187:           dir->col =0;
188:         }
189:         MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
190:         if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
191:           ISDestroy(dir->col);
192:           dir->col=0;
193:         }
194:         PetscOptionsHasName(pc->prefix,"-pc_cholesky_nonzeros_along_diagonal",&flg);
195:         if (flg) {
196:           PetscReal tol = 1.e-10;
197:           PetscOptionsGetReal(pc->prefix,"-pc_cholesky_nonzeros_along_diagonal",&tol,PETSC_NULL);
198:           MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
199:         }
200:         if (dir->row) {PetscLogObjectParent(pc,dir->row);}
201:       }
202:       MatDestroy(dir->fact);
203:       MatCholeskyFactorSymbolic(pc->pmat,dir->row,dir->info.fill,&dir->fact);
204:       MatGetInfo(dir->fact,MAT_LOCAL,&info);
205:       dir->actualfill = info.fill_ratio_needed;
206:       PetscLogObjectParent(pc,dir->fact);
207:     }
208:     MatCholeskyFactorNumeric(pc->pmat,&dir->fact);
209:   }
210:   return(0);
211: }

213: #undef __FUNCT__  
215: static int PCDestroy_Cholesky(PC pc)
216: {
217:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;
218:   int   ierr;

221:   if (!dir->inplace && dir->fact) {MatDestroy(dir->fact);}
222:   if (dir->row) {ISDestroy(dir->row);}
223:   if (dir->col) {ISDestroy(dir->col);}
224:   PetscStrfree(dir->ordering);
225:   PetscFree(dir);
226:   return(0);
227: }

229: #undef __FUNCT__  
231: static int PCApply_Cholesky(PC pc,Vec x,Vec y)
232: {
233:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;
234:   int   ierr;
235: 
237:   if (dir->inplace) {MatSolve(pc->pmat,x,y);}
238:   else              {MatSolve(dir->fact,x,y);}
239:   return(0);
240: }

242: #undef __FUNCT__  
244: static int PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
245: {
246:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;
247:   int   ierr;

250:   if (dir->inplace) {MatSolveTranspose(pc->pmat,x,y);}
251:   else              {MatSolveTranspose(dir->fact,x,y);}
252:   return(0);
253: }

255: /* -----------------------------------------------------------------------------------*/

257: EXTERN_C_BEGIN
258: #undef __FUNCT__  
260: int PCCholeskySetFill_Cholesky(PC pc,PetscReal fill)
261: {
262:   PC_Cholesky *dir;
263: 
265:   dir = (PC_Cholesky*)pc->data;
266:   dir->info.fill = fill;
267:   return(0);
268: }
269: EXTERN_C_END

271: EXTERN_C_BEGIN
272: #undef __FUNCT__  
274: int PCCholeskySetDamping_Cholesky(PC pc,PetscReal damping)
275: {
276:   PC_Cholesky *dir;
277: 
279:   dir = (PC_Cholesky*)pc->data;
280:   if (damping == (PetscReal) PETSC_DECIDE) {
281:     dir->info.damping = 1.e-12;
282:   } else {
283:     dir->info.damping = damping;
284:   }
285:   return(0);
286: }
287: EXTERN_C_END

289: EXTERN_C_BEGIN
290: #undef __FUNCT__  
292: int PCCholeskySetUseInPlace_Cholesky(PC pc)
293: {
294:   PC_Cholesky *dir;

297:   dir = (PC_Cholesky*)pc->data;
298:   dir->inplace = 1;
299:   return(0);
300: }
301: EXTERN_C_END

303: EXTERN_C_BEGIN
304: #undef __FUNCT__  
306: int PCCholeskySetMatOrdering_Cholesky(PC pc,MatOrderingType ordering)
307: {
308:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;
309:   int   ierr;

312:   PetscStrfree(dir->ordering);
313:   PetscStrallocpy(ordering,&dir->ordering);
314:   return(0);
315: }
316: EXTERN_C_END

318: /* -----------------------------------------------------------------------------------*/

320: #undef __FUNCT__  
322: /*@
323:    PCCholeskySetReuseOrdering - When similar matrices are factored, this
324:    causes the ordering computed in the first factor to be used for all
325:    following factors.

327:    Collective on PC

329:    Input Parameters:
330: +  pc - the preconditioner context
331: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

333:    Options Database Key:
334: .  -pc_cholesky_reuse_ordering - Activate PCCholeskySetReuseOrdering()

336:    Level: intermediate

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

340: .seealso: PCCholeskySetReuseFill(), PCICholeskySetReuseOrdering(), PCICholeskyDTSetReuseFill()
341: @*/
342: int PCCholeskySetReuseOrdering(PC pc,PetscTruth flag)
343: {
344:   int ierr,(*f)(PC,PetscTruth);

348:   PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetReuseOrdering_C",(void (**)(void))&f);
349:   if (f) {
350:     (*f)(pc,flag);
351:   }
352:   return(0);
353: }

355: #undef __FUNCT__  
357: /*@
358:    PCCholeskySetReuseFill - When matrices with same nonzero structure are Cholesky factored,
359:    this causes later ones to use the fill computed in the initial factorization.

361:    Collective on PC

363:    Input Parameters:
364: +  pc - the preconditioner context
365: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

367:    Options Database Key:
368: .  -pc_cholesky_reuse_fill - Activates PCCholeskySetReuseFill()

370:    Level: intermediate

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

374: .seealso: PCICholeskySetReuseOrdering(), PCCholeskySetReuseOrdering(), PCICholeskyDTSetReuseFill()
375: @*/
376: int PCCholeskySetReuseFill(PC pc,PetscTruth flag)
377: {
378:   int ierr,(*f)(PC,PetscTruth);

382:   PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetReuseFill_C",(void (**)(void))&f);
383:   if (f) {
384:     (*f)(pc,flag);
385:   }
386:   return(0);
387: }

389: #undef __FUNCT__  
391: /*@
392:    PCCholeskySetFill - Indicates the amount of fill you expect in the factored matrix,
393:    fill = number nonzeros in factor/number nonzeros in original matrix.

395:    Collective on PC
396:    
397:    Input Parameters:
398: +  pc - the preconditioner context
399: -  fill - amount of expected fill

401:    Options Database Key:
402: .  -pc_cholesky_fill <fill> - Sets fill amount

404:    Level: intermediate

406:    Note:
407:    For sparse matrix factorizations it is difficult to predict how much 
408:    fill to expect. By running with the option -log_info PETSc will print the 
409:    actual amount of fill used; allowing you to set the value accurately for
410:    future runs. Default PETSc uses a value of 5.0

412: .keywords: PC, set, factorization, direct, fill

414: .seealso: PCILUSetFill()
415: @*/
416: int PCCholeskySetFill(PC pc,PetscReal fill)
417: {
418:   int ierr,(*f)(PC,PetscReal);

422:   if (fill < 1.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Fill factor cannot be less then 1.0");
423:   PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetFill_C",(void (**)(void))&f);
424:   if (f) {
425:     (*f)(pc,fill);
426:   }
427:   return(0);
428: }

430: #undef __FUNCT__  
432: /*@
433:    PCCholeskySetDamping - Adds this quantity to the diagonal of the matrix during the 
434:    Cholesky numerical factorization.

436:    Collective on PC
437:    
438:    Input Parameters:
439: +  pc - the preconditioner context
440: -  damping - amount of damping

442:    Options Database Key:
443: .  -pc_cholesky_damping <damping> - Sets damping amount

445:    Level: intermediate

447: .keywords: PC, set, factorization, direct, fill

449: .seealso: PCICholeskySetFill(), PCILUSetDamp()
450: @*/
451: int PCCholeskySetDamping(PC pc,PetscReal damping)
452: {
453:   int ierr,(*f)(PC,PetscReal);

457:   PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetDamping_C",(void (**)(void))&f);
458:   if (f) {
459:     (*f)(pc,damping);
460:   }
461:   return(0);
462: }

464: #undef __FUNCT__  
466: /*@
467:    PCCholeskySetUseInPlace - Tells the system to do an in-place factorization.
468:    For dense matrices, this enables the solution of much larger problems. 
469:    For sparse matrices the factorization cannot be done truly in-place 
470:    so this does not save memory during the factorization, but after the matrix
471:    is factored, the original unfactored matrix is freed, thus recovering that
472:    space.

474:    Collective on PC

476:    Input Parameters:
477: .  pc - the preconditioner context

479:    Options Database Key:
480: .  -pc_cholesky_in_place - Activates in-place factorization

482:    Notes:
483:    PCCholeskySetUseInplace() can only be used with the KSP method KSPPREONLY or when 
484:    a different matrix is provided for the multiply and the preconditioner in 
485:    a call to SLESSetOperators().
486:    This is because the Krylov space methods require an application of the 
487:    matrix multiplication, which is not possible here because the matrix has 
488:    been factored in-place, replacing the original matrix.

490:    Level: intermediate

492: .keywords: PC, set, factorization, direct, inplace, in-place, Cholesky

494: .seealso: PCICholeskySetUseInPlace()
495: @*/
496: int PCCholeskySetUseInPlace(PC pc)
497: {
498:   int ierr,(*f)(PC);

502:   PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetUseInPlace_C",(void (**)(void))&f);
503:   if (f) {
504:     (*f)(pc);
505:   }
506:   return(0);
507: }

509: #undef __FUNCT__  
511: /*@
512:     PCCholeskySetMatOrdering - Sets the ordering routine (to reduce fill) to 
513:     be used it the Cholesky factorization.

515:     Collective on PC

517:     Input Parameters:
518: +   pc - the preconditioner context
519: -   ordering - the matrix ordering name, for example, MATORDERING_ND or MATORDERING_RCM

521:     Options Database Key:
522: .   -pc_cholesky_mat_ordering_type <nd,rcm,...> - Sets ordering routine

524:     Level: intermediate

526: .seealso: PCICholeskySetMatOrdering()
527: @*/
528: int PCCholeskySetMatOrdering(PC pc,MatOrderingType ordering)
529: {
530:   int ierr,(*f)(PC,MatOrderingType);

533:   PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetMatOrdering_C",(void (**)(void))&f);
534:   if (f) {
535:     (*f)(pc,ordering);
536:   }
537:   return(0);
538: }

540: /*MC
541:    PCCholesky - Uses a direct solver, based on Cholesky factorization, as a preconditioner

543:    Options Database Keys:
544: +  -pc_cholesky_reuse_ordering - Activate PCLUSetReuseOrdering()
545: .  -pc_cholesky_reuse_fill - Activates PCLUSetReuseFill()
546: .  -pc_cholesky_fill <fill> - Sets fill amount
547: .  -pc_cholesky_damping <damping> - Sets damping amount
548: .  -pc_cholesky_in_place - Activates in-place factorization
549: -  -pc_cholesky_mat_ordering_type <nd,rcm,...> - Sets ordering routine

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

553:    Level: beginner

555:    Concepts: Cholesky factorization, direct solver

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

561: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
562:            PCILU, PCLU, PCICC, PCCholeskySetReuseOrdering(), PCCholeskySetReuseFill(), PCGetFactoredMatrix(),
563:            PCCholeskySetFill(), PCCholeskySetDamping(), PCCholeskySetUseInPlace(), PCCholeskySetMatOrdering()

565: M*/

567: EXTERN_C_BEGIN
568: #undef __FUNCT__  
570: int PCCreate_Cholesky(PC pc)
571: {
572:   int         ierr;
573:   PC_Cholesky *dir;

576:   PetscNew(PC_Cholesky,&dir);
577:   PetscLogObjectMemory(pc,sizeof(PC_Cholesky));

579:   dir->fact               = 0;
580:   dir->inplace            = 0;
581:   dir->info.fill          = 5.0;
582:   dir->info.damping       = 0.0;
583:   dir->info.pivotinblocks = 1.0;
584:   dir->col                = 0;
585:   dir->row                = 0;
586:   PetscStrallocpy(MATORDERING_NATURAL,&dir->ordering);
587:   dir->reusefill        = PETSC_FALSE;
588:   dir->reuseordering    = PETSC_FALSE;
589:   pc->data              = (void*)dir;

591:   pc->ops->destroy           = PCDestroy_Cholesky;
592:   pc->ops->apply             = PCApply_Cholesky;
593:   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
594:   pc->ops->setup             = PCSetUp_Cholesky;
595:   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
596:   pc->ops->view              = PCView_Cholesky;
597:   pc->ops->applyrichardson   = 0;
598:   pc->ops->getfactoredmatrix = PCGetFactoredMatrix_Cholesky;

600:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetFill_C","PCCholeskySetFill_Cholesky",
601:                     PCCholeskySetFill_Cholesky);
602:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetDamping_C","PCCholeskySetDamping_Cholesky",
603:                     PCCholeskySetDamping_Cholesky);
604:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetUseInPlace_C","PCCholeskySetUseInPlace_Cholesky",
605:                     PCCholeskySetUseInPlace_Cholesky);
606:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetMatOrdering_C","PCCholeskySetMatOrdering_Cholesky",
607:                     PCCholeskySetMatOrdering_Cholesky);
608:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetReuseOrdering_C","PCCholeskySetReuseOrdering_Cholesky",
609:                     PCCholeskySetReuseOrdering_Cholesky);
610:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetReuseFill_C","PCCholeskySetReuseFill_Cholesky",
611:                     PCCholeskySetReuseFill_Cholesky);
612:   return(0);
613: }
614: EXTERN_C_END