Actual source code: eisen.c
2: /*
3: Defines a Eisenstat trick SSOR preconditioner. This uses about
4: %50 of the usual amount of floating point ops used for SSOR + Krylov
5: method. But it requires actually solving the preconditioned problem
6: with both left and right preconditioning.
7: */
8: #include <petsc/private/pcimpl.h>
10: typedef struct {
11: Mat shell,A;
12: Vec b[2],diag; /* temporary storage for true right hand side */
13: PetscReal omega;
14: PetscBool usediag; /* indicates preconditioner should include diagonal scaling*/
15: } PC_Eisenstat;
17: static PetscErrorCode PCMult_Eisenstat(Mat mat,Vec b,Vec x)
18: {
19: PC pc;
20: PC_Eisenstat *eis;
22: MatShellGetContext(mat,&pc);
23: eis = (PC_Eisenstat*)pc->data;
24: MatSOR(eis->A,b,eis->omega,SOR_EISENSTAT,0.0,1,1,x);
25: return 0;
26: }
28: static PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y)
29: {
30: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
31: PetscBool hasop;
33: if (eis->usediag) {
34: MatHasOperation(pc->pmat,MATOP_MULT_DIAGONAL_BLOCK,&hasop);
35: if (hasop) {
36: MatMultDiagonalBlock(pc->pmat,x,y);
37: } else {
38: VecPointwiseMult(y,x,eis->diag);
39: }
40: } else VecCopy(x,y);
41: return 0;
42: }
44: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
45: {
46: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
47: PetscBool nonzero;
49: if (pc->presolvedone < 2) {
51: /* swap shell matrix and true matrix */
52: eis->A = pc->mat;
53: pc->mat = eis->shell;
54: }
56: if (!eis->b[pc->presolvedone-1]) {
57: VecDuplicate(b,&eis->b[pc->presolvedone-1]);
58: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->b[pc->presolvedone-1]);
59: }
61: /* if nonzero initial guess, modify x */
62: KSPGetInitialGuessNonzero(ksp,&nonzero);
63: if (nonzero) {
64: VecCopy(x,eis->b[pc->presolvedone-1]);
65: MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
66: }
68: /* save true b, other option is to swap pointers */
69: VecCopy(b,eis->b[pc->presolvedone-1]);
71: /* modify b by (L + D/omega)^{-1} */
72: MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP),0.0,1,1,b);
73: return 0;
74: }
76: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
77: {
78: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
80: /* get back true b */
81: VecCopy(eis->b[pc->presolvedone],b);
83: /* modify x by (U + D/omega)^{-1} */
84: VecCopy(x,eis->b[pc->presolvedone]);
85: MatSOR(eis->A,eis->b[pc->presolvedone],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP),0.0,1,1,x);
86: if (!pc->presolvedone) pc->mat = eis->A;
87: return 0;
88: }
90: static PetscErrorCode PCReset_Eisenstat(PC pc)
91: {
92: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
94: VecDestroy(&eis->b[0]);
95: VecDestroy(&eis->b[1]);
96: MatDestroy(&eis->shell);
97: VecDestroy(&eis->diag);
98: return 0;
99: }
101: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
102: {
103: PCReset_Eisenstat(pc);
104: PetscFree(pc->data);
105: return 0;
106: }
108: static PetscErrorCode PCSetFromOptions_Eisenstat(PetscOptionItems *PetscOptionsObject,PC pc)
109: {
110: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
111: PetscBool set,flg;
113: PetscOptionsHead(PetscOptionsObject,"Eisenstat SSOR options");
114: PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,NULL);
115: PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatSetNoDiagonalScaling",eis->usediag ? PETSC_FALSE : PETSC_TRUE,&flg,&set);
116: if (set) {
117: PCEisenstatSetNoDiagonalScaling(pc,flg);
118: }
119: PetscOptionsTail();
120: return 0;
121: }
123: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
124: {
125: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
126: PetscBool iascii;
128: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
129: if (iascii) {
130: PetscViewerASCIIPrintf(viewer," omega = %g\n",(double)eis->omega);
131: if (eis->usediag) {
132: PetscViewerASCIIPrintf(viewer," Using diagonal scaling (default)\n");
133: } else {
134: PetscViewerASCIIPrintf(viewer," Not using diagonal scaling\n");
135: }
136: }
137: return 0;
138: }
140: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
141: {
142: PetscInt M,N,m,n;
143: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
145: if (!pc->setupcalled) {
146: MatGetSize(pc->mat,&M,&N);
147: MatGetLocalSize(pc->mat,&m,&n);
148: MatCreate(PetscObjectComm((PetscObject)pc),&eis->shell);
149: MatSetSizes(eis->shell,m,n,M,N);
150: MatSetType(eis->shell,MATSHELL);
151: MatSetUp(eis->shell);
152: MatShellSetContext(eis->shell,pc);
153: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->shell);
154: MatShellSetOperation(eis->shell,MATOP_MULT,(void (*)(void))PCMult_Eisenstat);
155: }
156: if (!eis->usediag) return 0;
157: if (!pc->setupcalled) {
158: MatCreateVecs(pc->pmat,&eis->diag,NULL);
159: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->diag);
160: }
161: MatGetDiagonal(pc->pmat,eis->diag);
162: return 0;
163: }
165: /* --------------------------------------------------------------------*/
167: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
168: {
169: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
172: eis->omega = omega;
173: return 0;
174: }
176: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc,PetscBool flg)
177: {
178: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
180: eis->usediag = flg;
181: return 0;
182: }
184: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc,PetscReal *omega)
185: {
186: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
188: *omega = eis->omega;
189: return 0;
190: }
192: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc,PetscBool *flg)
193: {
194: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
196: *flg = eis->usediag;
197: return 0;
198: }
200: /*@
201: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
202: to use with Eisenstat's trick (where omega = 1.0 by default).
204: Logically Collective on PC
206: Input Parameters:
207: + pc - the preconditioner context
208: - omega - relaxation coefficient (0 < omega < 2)
210: Options Database Key:
211: . -pc_eisenstat_omega <omega> - Sets omega
213: Notes:
214: The Eisenstat trick implementation of SSOR requires about 50% of the
215: usual amount of floating point operations used for SSOR + Krylov method;
216: however, the preconditioned problem must be solved with both left
217: and right preconditioning.
219: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
220: which can be chosen with the database options
221: $ -pc_type sor -pc_sor_symmetric
223: Level: intermediate
225: .seealso: PCSORSetOmega()
226: @*/
227: PetscErrorCode PCEisenstatSetOmega(PC pc,PetscReal omega)
228: {
231: PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
232: return 0;
233: }
235: /*@
236: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner
237: not to do additional diagonal preconditioning. For matrices with a constant
238: along the diagonal, this may save a small amount of work.
240: Logically Collective on PC
242: Input Parameters:
243: + pc - the preconditioner context
244: - flg - PETSC_TRUE turns off diagonal scaling inside the algorithm
246: Options Database Key:
247: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
249: Level: intermediate
251: Note:
252: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
253: likley want to use this routine since it will save you some unneeded flops.
255: .seealso: PCEisenstatSetOmega()
256: @*/
257: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)
258: {
260: PetscTryMethod(pc,"PCEisenstatSetNoDiagonalScaling_C",(PC,PetscBool),(pc,flg));
261: return 0;
262: }
264: /*@
265: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
266: to use with Eisenstat's trick (where omega = 1.0 by default).
268: Logically Collective on PC
270: Input Parameter:
271: . pc - the preconditioner context
273: Output Parameter:
274: . omega - relaxation coefficient (0 < omega < 2)
276: Options Database Key:
277: . -pc_eisenstat_omega <omega> - Sets omega
279: Notes:
280: The Eisenstat trick implementation of SSOR requires about 50% of the
281: usual amount of floating point operations used for SSOR + Krylov method;
282: however, the preconditioned problem must be solved with both left
283: and right preconditioning.
285: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
286: which can be chosen with the database options
287: $ -pc_type sor -pc_sor_symmetric
289: Level: intermediate
291: .seealso: PCSORGetOmega(), PCEisenstatSetOmega()
292: @*/
293: PetscErrorCode PCEisenstatGetOmega(PC pc,PetscReal *omega)
294: {
296: PetscUseMethod(pc,"PCEisenstatGetOmega_C",(PC,PetscReal*),(pc,omega));
297: return 0;
298: }
300: /*@
301: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
302: not to do additional diagonal preconditioning. For matrices with a constant
303: along the diagonal, this may save a small amount of work.
305: Logically Collective on PC
307: Input Parameter:
308: . pc - the preconditioner context
310: Output Parameter:
311: . flg - PETSC_TRUE means there is no diagonal scaling applied
313: Options Database Key:
314: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
316: Level: intermediate
318: Note:
319: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
320: likley want to use this routine since it will save you some unneeded flops.
322: .seealso: PCEisenstatGetOmega()
323: @*/
324: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool *flg)
325: {
327: PetscUseMethod(pc,"PCEisenstatGetNoDiagonalScaling_C",(PC,PetscBool*),(pc,flg));
328: return 0;
329: }
331: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool* change)
332: {
333: *change = PETSC_TRUE;
334: return 0;
335: }
337: /* --------------------------------------------------------------------*/
339: /*MC
340: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
341: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
343: Options Database Keys:
344: + -pc_eisenstat_omega <omega> - Sets omega
345: - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
347: Level: beginner
349: Notes:
350: Only implemented for the SeqAIJ matrix format.
351: Not a true parallel SOR, in parallel this implementation corresponds to block
352: Jacobi with SOR on each block.
354: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
355: PCEisenstatSetNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
356: M*/
358: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
359: {
360: PC_Eisenstat *eis;
362: PetscNewLog(pc,&eis);
364: pc->ops->apply = PCApply_Eisenstat;
365: pc->ops->presolve = PCPreSolve_Eisenstat;
366: pc->ops->postsolve = PCPostSolve_Eisenstat;
367: pc->ops->applyrichardson = NULL;
368: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
369: pc->ops->destroy = PCDestroy_Eisenstat;
370: pc->ops->reset = PCReset_Eisenstat;
371: pc->ops->view = PCView_Eisenstat;
372: pc->ops->setup = PCSetUp_Eisenstat;
374: pc->data = eis;
375: eis->omega = 1.0;
376: eis->b[0] = NULL;
377: eis->b[1] = NULL;
378: eis->diag = NULL;
379: eis->usediag = PETSC_TRUE;
381: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetOmega_C",PCEisenstatSetOmega_Eisenstat);
382: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetNoDiagonalScaling_C",PCEisenstatSetNoDiagonalScaling_Eisenstat);
383: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetOmega_C",PCEisenstatGetOmega_Eisenstat);
384: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetNoDiagonalScaling_C",PCEisenstatGetNoDiagonalScaling_Eisenstat);
385: PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_Eisenstat);
386: return 0;
387: }