xref: /petsc/src/ksp/pc/impls/factor/cholesky/cholesky.c (revision b271bb04eaec0d477b629c8be871cf42cb4980f5)
1 #define PETSCKSP_DLL
2 
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"                /*I "petscpc.h" I*/
9 
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   MatSolverPackage solvertype;
20 } PC_Cholesky;
21 
22 EXTERN_C_BEGIN
23 #undef __FUNCT__
24 #define __FUNCT__ "PCFactorSetMatSolverPackage_Cholesky"
25 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatSolverPackage_Cholesky(PC pc,const MatSolverPackage stype)
26 {
27   PetscErrorCode ierr;
28   PC_Cholesky    *choleksy = (PC_Cholesky*)pc->data;
29 
30   PetscFunctionBegin;
31   ierr = PetscStrfree(choleksy->solvertype);CHKERRQ(ierr);
32   ierr = PetscStrallocpy(stype,&choleksy->solvertype);CHKERRQ(ierr);
33   PetscFunctionReturn(0);
34 }
35 EXTERN_C_END
36 
37 EXTERN_C_BEGIN
38 #undef __FUNCT__
39 #define __FUNCT__ "PCFactorSetZeroPivot_Cholesky"
40 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetZeroPivot_Cholesky(PC pc,PetscReal z)
41 {
42   PC_Cholesky *ch;
43 
44   PetscFunctionBegin;
45   ch                 = (PC_Cholesky*)pc->data;
46   ch->info.zeropivot = z;
47   PetscFunctionReturn(0);
48 }
49 EXTERN_C_END
50 
51 EXTERN_C_BEGIN
52 #undef __FUNCT__
53 #define __FUNCT__ "PCFactorSetShiftNonzero_Cholesky"
54 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftNonzero_Cholesky(PC pc,PetscReal shift)
55 {
56   PC_Cholesky *dir;
57 
58   PetscFunctionBegin;
59   dir = (PC_Cholesky*)pc->data;
60   if (shift == (PetscReal) PETSC_DECIDE) {
61     dir->info.shiftnz = 1.e-12;
62   } else {
63     dir->info.shiftnz = shift;
64   }
65   PetscFunctionReturn(0);
66 }
67 EXTERN_C_END
68 
69 EXTERN_C_BEGIN
70 #undef __FUNCT__
71 #define __FUNCT__ "PCFactorSetShiftPd_Cholesky"
72 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftPd_Cholesky(PC pc,PetscTruth shift)
73 {
74   PC_Cholesky *dir;
75 
76   PetscFunctionBegin;
77   dir = (PC_Cholesky*)pc->data;
78   if (shift) {
79     dir->info.shiftpd = 1.0;
80   } else {
81     dir->info.shiftpd = 0.0;
82   }
83   PetscFunctionReturn(0);
84 }
85 EXTERN_C_END
86 
87 EXTERN_C_BEGIN
88 #undef __FUNCT__
89 #define __FUNCT__ "PCFactorSetReuseOrdering_Cholesky"
90 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering_Cholesky(PC pc,PetscTruth flag)
91 {
92   PC_Cholesky *lu;
93 
94   PetscFunctionBegin;
95   lu               = (PC_Cholesky*)pc->data;
96   lu->reuseordering = flag;
97   PetscFunctionReturn(0);
98 }
99 EXTERN_C_END
100 
101 EXTERN_C_BEGIN
102 #undef __FUNCT__
103 #define __FUNCT__ "PCFactorSetReuseFill_Cholesky"
104 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseFill_Cholesky(PC pc,PetscTruth flag)
105 {
106   PC_Cholesky *lu;
107 
108   PetscFunctionBegin;
109   lu = (PC_Cholesky*)pc->data;
110   lu->reusefill = flag;
111   PetscFunctionReturn(0);
112 }
113 EXTERN_C_END
114 
115 #undef __FUNCT__
116 #define __FUNCT__ "PCSetFromOptions_Cholesky"
117 static PetscErrorCode PCSetFromOptions_Cholesky(PC pc)
118 {
119   PC_Cholesky    *lu = (PC_Cholesky*)pc->data;
120   PetscErrorCode ierr;
121   PetscTruth     flg;
122   char           tname[256], solvertype[64];
123   PetscFList     ordlist;
124 
125   PetscFunctionBegin;
126   ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr);
127   ierr = PetscOptionsHead("Cholesky options");CHKERRQ(ierr);
128     ierr = PetscOptionsName("-pc_factor_in_place","Form Cholesky in the same memory as the matrix","PCFactorSetUseInPlace",&flg);CHKERRQ(ierr);
129     if (flg) {
130       ierr = PCFactorSetUseInPlace(pc);CHKERRQ(ierr);
131     }
132     ierr = PetscOptionsReal("-pc_factor_fill","Expected non-zeros in Cholesky/non-zeros in matrix","PCFactorSetFill",lu->info.fill,&lu->info.fill,0);CHKERRQ(ierr);
133 
134     ierr = PetscOptionsName("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",&flg);CHKERRQ(ierr);
135     if (flg) {
136       ierr = PCFactorSetReuseFill(pc,PETSC_TRUE);CHKERRQ(ierr);
137     }
138     ierr = PetscOptionsName("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",&flg);CHKERRQ(ierr);
139     if (flg) {
140       ierr = PCFactorSetReuseOrdering(pc,PETSC_TRUE);CHKERRQ(ierr);
141     }
142 
143     ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr);
144     ierr = PetscOptionsList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in Cholesky","PCFactorSetMatOrderingType",ordlist,lu->ordering,tname,256,&flg);CHKERRQ(ierr);
145     if (flg) {
146       ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr);
147     }
148 
149     /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
150     ierr = PetscOptionsString("-pc_factor_mat_solver_package","Specific Cholesky solver to use","MatGetFactor",lu->solvertype,solvertype,64,&flg);CHKERRQ(ierr);
151     if (flg) {
152       ierr = PCFactorSetMatSolverPackage(pc,solvertype);CHKERRQ(ierr);
153     }
154 
155     ierr = PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);CHKERRQ(ierr);
156     if (flg) {
157       ierr = PCFactorSetShiftNonzero(pc,(PetscReal) PETSC_DECIDE);CHKERRQ(ierr);
158     }
159     ierr = PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",lu->info.shiftnz,&lu->info.shiftnz,0);CHKERRQ(ierr);
160     ierr = PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",&flg);CHKERRQ(ierr);
161     if (flg) {
162       ierr = PCFactorSetShiftPd(pc,PETSC_TRUE);CHKERRQ(ierr);
163     }
164     ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",lu->info.zeropivot,&lu->info.zeropivot,0);CHKERRQ(ierr);
165 
166   ierr = PetscOptionsTail();CHKERRQ(ierr);
167   PetscFunctionReturn(0);
168 }
169 
170 #undef __FUNCT__
171 #define __FUNCT__ "PCView_Cholesky"
172 static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer)
173 {
174   PC_Cholesky    *lu = (PC_Cholesky*)pc->data;
175   PetscErrorCode ierr;
176   PetscTruth     iascii,isstring;
177 
178   PetscFunctionBegin;
179   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
180   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr);
181   if (iascii) {
182 
183     if (lu->inplace) {ierr = PetscViewerASCIIPrintf(viewer,"  Cholesky: in-place factorization\n");CHKERRQ(ierr);}
184     else             {ierr = PetscViewerASCIIPrintf(viewer,"  Cholesky: out-of-place factorization\n");CHKERRQ(ierr);}
185     ierr = PetscViewerASCIIPrintf(viewer,"    matrix ordering: %s\n",lu->ordering);CHKERRQ(ierr);
186     if (lu->fact) {
187       ierr = PetscViewerASCIIPrintf(viewer,"  Cholesky: factor fill ratio needed %G\n",lu->actualfill);CHKERRQ(ierr);
188       ierr = PetscViewerASCIIPrintf(viewer,"       Factored matrix follows\n");CHKERRQ(ierr);
189       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
190       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
191       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
192       ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
193       ierr = MatView(lu->fact,viewer);CHKERRQ(ierr);
194       ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
195       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
196       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
197       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
198     }
199     if (lu->reusefill)    {ierr = PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorization\n");CHKERRQ(ierr);}
200     if (lu->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorization\n");CHKERRQ(ierr);}
201   } else if (isstring) {
202     ierr = PetscViewerStringSPrintf(viewer," order=%s",lu->ordering);CHKERRQ(ierr);CHKERRQ(ierr);
203   } else {
204     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCCHOLESKY",((PetscObject)viewer)->type_name);
205   }
206   PetscFunctionReturn(0);
207 }
208 
209 #undef __FUNCT__
210 #define __FUNCT__ "PCFactorGetMatrix_Cholesky"
211 static PetscErrorCode PCFactorGetMatrix_Cholesky(PC pc,Mat *mat)
212 {
213   PC_Cholesky *dir = (PC_Cholesky*)pc->data;
214 
215   PetscFunctionBegin;
216   if (!dir->fact) SETERRQ(PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
217   *mat = dir->fact;
218   PetscFunctionReturn(0);
219 }
220 
221 #undef __FUNCT__
222 #define __FUNCT__ "PCSetUp_Cholesky"
223 static PetscErrorCode PCSetUp_Cholesky(PC pc)
224 {
225   PetscErrorCode ierr;
226   PetscTruth     flg;
227   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
228 
229   PetscFunctionBegin;
230   if (dir->reusefill && pc->setupcalled) dir->info.fill = dir->actualfill;
231 
232   if (dir->inplace) {
233     if (dir->row && dir->col && (dir->row != dir->col)) {
234       ierr = ISDestroy(dir->row);CHKERRQ(ierr);
235       dir->row = 0;
236     }
237     if (dir->col) {
238       ierr = ISDestroy(dir->col);CHKERRQ(ierr);
239       dir->col = 0;
240     }
241     ierr = MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);CHKERRQ(ierr);
242     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
243       ierr = ISDestroy(dir->col);CHKERRQ(ierr);
244       dir->col=0;
245     }
246     if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);}
247     ierr = MatCholeskyFactor(pc->pmat,dir->row,&dir->info);CHKERRQ(ierr);
248     dir->fact = pc->pmat;
249   } else {
250     MatInfo info;
251     if (!pc->setupcalled) {
252       ierr = MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);CHKERRQ(ierr);
253       /* check if dir->row == dir->col */
254       ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr);
255       if (!flg) SETERRQ(PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
256       ierr = ISDestroy(dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */
257       dir->col=0;
258 
259       ierr = PetscOptionsHasName(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);CHKERRQ(ierr);
260       if (flg) {
261         PetscReal tol = 1.e-10;
262         ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr);
263         ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr);
264       }
265       if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);}
266       ierr = MatGetFactor(pc->pmat,dir->solvertype,MAT_FACTOR_CHOLESKY,&dir->fact);CHKERRQ(ierr);
267       ierr = MatCholeskyFactorSymbolic(dir->fact,pc->pmat,dir->row,&dir->info);CHKERRQ(ierr);
268       ierr = MatGetInfo(dir->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
269       dir->actualfill = info.fill_ratio_needed;
270       ierr = PetscLogObjectParent(pc,dir->fact);CHKERRQ(ierr);
271     } else if (pc->flag != SAME_NONZERO_PATTERN) {
272       if (!dir->reuseordering) {
273         if (dir->row && dir->col && (dir->row != dir->col)) {
274           ierr = ISDestroy(dir->row);CHKERRQ(ierr);
275           dir->row = 0;
276         }
277         if (dir->col) {
278           ierr = ISDestroy(dir->col);CHKERRQ(ierr);
279           dir->col =0;
280         }
281         ierr = MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);CHKERRQ(ierr);
282         if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
283           ierr = ISDestroy(dir->col);CHKERRQ(ierr);
284           dir->col=0;
285         }
286         ierr = PetscOptionsHasName(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);CHKERRQ(ierr);
287         if (flg) {
288           PetscReal tol = 1.e-10;
289           ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr);
290           ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr);
291         }
292         if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);}
293       }
294       ierr = MatDestroy(dir->fact);CHKERRQ(ierr);
295       ierr = MatGetFactor(pc->pmat,dir->solvertype,MAT_FACTOR_CHOLESKY,&dir->fact);CHKERRQ(ierr);
296       ierr = MatCholeskyFactorSymbolic(dir->fact,pc->pmat,dir->row,&dir->info);CHKERRQ(ierr);
297       ierr = MatGetInfo(dir->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
298       dir->actualfill = info.fill_ratio_needed;
299       ierr = PetscLogObjectParent(pc,dir->fact);CHKERRQ(ierr);
300     }
301     ierr = MatCholeskyFactorNumeric(dir->fact,pc->pmat,&dir->info);CHKERRQ(ierr);
302   }
303   PetscFunctionReturn(0);
304 }
305 
306 #undef __FUNCT__
307 #define __FUNCT__ "PCDestroy_Cholesky"
308 static PetscErrorCode PCDestroy_Cholesky(PC pc)
309 {
310   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
311   PetscErrorCode ierr;
312 
313   PetscFunctionBegin;
314   if (!dir->inplace && dir->fact) {ierr = MatDestroy(dir->fact);CHKERRQ(ierr);}
315   if (dir->row) {ierr = ISDestroy(dir->row);CHKERRQ(ierr);}
316   if (dir->col) {ierr = ISDestroy(dir->col);CHKERRQ(ierr);}
317   ierr = PetscStrfree(dir->ordering);CHKERRQ(ierr);
318   ierr = PetscStrfree(dir->solvertype);CHKERRQ(ierr);
319   ierr = PetscFree(dir);CHKERRQ(ierr);
320   PetscFunctionReturn(0);
321 }
322 
323 #undef __FUNCT__
324 #define __FUNCT__ "PCApply_Cholesky"
325 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
326 {
327   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
328   PetscErrorCode ierr;
329 
330   PetscFunctionBegin;
331   if (dir->inplace) {ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr);}
332   else              {ierr = MatSolve(dir->fact,x,y);CHKERRQ(ierr);}
333   PetscFunctionReturn(0);
334 }
335 
336 #undef __FUNCT__
337 #define __FUNCT__ "PCApplyTranspose_Cholesky"
338 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
339 {
340   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
341   PetscErrorCode ierr;
342 
343   PetscFunctionBegin;
344   if (dir->inplace) {ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr);}
345   else              {ierr = MatSolveTranspose(dir->fact,x,y);CHKERRQ(ierr);}
346   PetscFunctionReturn(0);
347 }
348 
349 /* -----------------------------------------------------------------------------------*/
350 
351 EXTERN_C_BEGIN
352 #undef __FUNCT__
353 #define __FUNCT__ "PCFactorSetFill_Cholesky"
354 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetFill_Cholesky(PC pc,PetscReal fill)
355 {
356   PC_Cholesky *dir;
357 
358   PetscFunctionBegin;
359   dir = (PC_Cholesky*)pc->data;
360   dir->info.fill = fill;
361   PetscFunctionReturn(0);
362 }
363 EXTERN_C_END
364 
365 EXTERN_C_BEGIN
366 #undef __FUNCT__
367 #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky"
368 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseInPlace_Cholesky(PC pc)
369 {
370   PC_Cholesky *dir;
371 
372   PetscFunctionBegin;
373   dir = (PC_Cholesky*)pc->data;
374   dir->inplace = PETSC_TRUE;
375   PetscFunctionReturn(0);
376 }
377 EXTERN_C_END
378 
379 EXTERN_C_BEGIN
380 #undef __FUNCT__
381 #define __FUNCT__ "PCFactorSetMatOrderingType_Cholesky"
382 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatOrderingType_Cholesky(PC pc,const MatOrderingType ordering)
383 {
384   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
385   PetscErrorCode ierr;
386 
387   PetscFunctionBegin;
388   ierr = PetscStrfree(dir->ordering);CHKERRQ(ierr);
389   ierr = PetscStrallocpy(ordering,&dir->ordering);CHKERRQ(ierr);
390   PetscFunctionReturn(0);
391 }
392 EXTERN_C_END
393 
394 /* -----------------------------------------------------------------------------------*/
395 
396 #undef __FUNCT__
397 #define __FUNCT__ "PCFactorSetReuseOrdering"
398 /*@
399    PCFactorSetReuseOrdering - When similar matrices are factored, this
400    causes the ordering computed in the first factor to be used for all
401    following factors.
402 
403    Collective on PC
404 
405    Input Parameters:
406 +  pc - the preconditioner context
407 -  flag - PETSC_TRUE to reuse else PETSC_FALSE
408 
409    Options Database Key:
410 .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
411 
412    Level: intermediate
413 
414 .keywords: PC, levels, reordering, factorization, incomplete, LU
415 
416 .seealso: PCFactorSetReuseFill()
417 @*/
418 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering(PC pc,PetscTruth flag)
419 {
420   PetscErrorCode ierr,(*f)(PC,PetscTruth);
421 
422   PetscFunctionBegin;
423   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
424   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",(void (**)(void))&f);CHKERRQ(ierr);
425   if (f) {
426     ierr = (*f)(pc,flag);CHKERRQ(ierr);
427   }
428   PetscFunctionReturn(0);
429 }
430 
431 #undef __FUNCT__
432 #define __FUNCT__ "PCFactorSetReuseFill"
433 /*@
434    PCFactorSetReuseFill - When matrices with same different nonzero structure are factored,
435    this causes later ones to use the fill ratio computed in the initial factorization.
436 
437    Collective on PC
438 
439    Input Parameters:
440 +  pc - the preconditioner context
441 -  flag - PETSC_TRUE to reuse else PETSC_FALSE
442 
443    Options Database Key:
444 .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
445 
446    Level: intermediate
447 
448 .keywords: PC, levels, reordering, factorization, incomplete, Cholesky
449 
450 .seealso: PCFactorSetReuseOrdering()
451 @*/
452 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseFill(PC pc,PetscTruth flag)
453 {
454   PetscErrorCode ierr,(*f)(PC,PetscTruth);
455 
456   PetscFunctionBegin;
457   PetscValidHeaderSpecific(pc,PC_COOKIE,2);
458   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseFill_C",(void (**)(void))&f);CHKERRQ(ierr);
459   if (f) {
460     ierr = (*f)(pc,flag);CHKERRQ(ierr);
461   }
462   PetscFunctionReturn(0);
463 }
464 
465 /*MC
466    PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
467 
468    Options Database Keys:
469 +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
470 .  -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like spooles
471 .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
472 .  -pc_factor_fill <fill> - Sets fill amount
473 .  -pc_factor_in_place - Activates in-place factorization
474 .  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
475 .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
476 -  -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
477    is optional with PETSC_TRUE being the default
478 
479    Notes: Not all options work for all matrix formats
480 
481    Level: beginner
482 
483    Concepts: Cholesky factorization, direct solver
484 
485    Notes: Usually this will compute an "exact" solution in one iteration and does
486           not need a Krylov method (i.e. you can use -ksp_type preonly, or
487           KSPSetType(ksp,KSPPREONLY) for the Krylov method
488 
489 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
490            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
491            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(),
492 	   PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd()
493 
494 M*/
495 
496 EXTERN_C_BEGIN
497 #undef __FUNCT__
498 #define __FUNCT__ "PCCreate_Cholesky"
499 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Cholesky(PC pc)
500 {
501   PetscErrorCode ierr;
502   PC_Cholesky    *dir;
503 
504   PetscFunctionBegin;
505   ierr = PetscNewLog(pc,PC_Cholesky,&dir);CHKERRQ(ierr);
506 
507   dir->fact                   = 0;
508   dir->inplace                = PETSC_FALSE;
509   ierr = MatFactorInfoInitialize(&dir->info);CHKERRQ(ierr);
510   dir->info.fill              = 5.0;
511   dir->info.shiftnz           = 0.0;
512   dir->info.shiftpd           = 0.0; /* false */
513   dir->info.pivotinblocks     = 1.0;
514   dir->col                    = 0;
515   dir->row                    = 0;
516   ierr = PetscStrallocpy(MATORDERING_NATURAL,&dir->ordering);CHKERRQ(ierr);
517   ierr = PetscStrallocpy(MAT_SOLVER_PETSC,&dir->solvertype);CHKERRQ(ierr);
518   dir->reusefill        = PETSC_FALSE;
519   dir->reuseordering    = PETSC_FALSE;
520   pc->data              = (void*)dir;
521 
522   pc->ops->destroy           = PCDestroy_Cholesky;
523   pc->ops->apply             = PCApply_Cholesky;
524   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
525   pc->ops->setup             = PCSetUp_Cholesky;
526   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
527   pc->ops->view              = PCView_Cholesky;
528   pc->ops->applyrichardson   = 0;
529   pc->ops->getfactoredmatrix = PCFactorGetMatrix_Cholesky;
530 
531   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Cholesky",
532                     PCFactorSetMatSolverPackage_Cholesky);CHKERRQ(ierr);
533   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Cholesky",
534                     PCFactorSetZeroPivot_Cholesky);CHKERRQ(ierr);
535   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Cholesky",
536                     PCFactorSetShiftNonzero_Cholesky);CHKERRQ(ierr);
537   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Cholesky",
538                     PCFactorSetShiftPd_Cholesky);CHKERRQ(ierr);
539 
540   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Cholesky",
541                     PCFactorSetFill_Cholesky);CHKERRQ(ierr);
542   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky",
543                     PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr);
544   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Cholesky",
545                     PCFactorSetMatOrderingType_Cholesky);CHKERRQ(ierr);
546   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky",
547                     PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr);
548   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky",
549                     PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr);
550   PetscFunctionReturn(0);
551 }
552 EXTERN_C_END
553