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: {
20: PC pc;
21: PC_Eisenstat *eis;
24: MatShellGetContext(mat,&pc);
25: eis = (PC_Eisenstat*)pc->data;
26: MatSOR(eis->A,b,eis->omega,SOR_EISENSTAT,0.0,1,1,x);
27: return(0);
28: }
30: static PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y)
31: {
32: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
34: PetscBool hasop;
37: if (eis->usediag) {
38: MatHasOperation(pc->pmat,MATOP_MULT_DIAGONAL_BLOCK,&hasop);
39: if (hasop) {
40: MatMultDiagonalBlock(pc->pmat,x,y);
41: } else {
42: VecPointwiseMult(y,x,eis->diag);
43: }
44: } else {VecCopy(x,y);}
45: return(0);
46: }
48: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
49: {
50: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
51: PetscBool nonzero;
55: if (pc->presolvedone < 2) {
56: if (pc->mat != pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot have different mat and pmat");
57: /* swap shell matrix and true matrix */
58: eis->A = pc->mat;
59: pc->mat = eis->shell;
60: }
62: if (!eis->b[pc->presolvedone-1]) {
63: VecDuplicate(b,&eis->b[pc->presolvedone-1]);
64: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->b[pc->presolvedone-1]);
65: }
67: /* if nonzero initial guess, modify x */
68: KSPGetInitialGuessNonzero(ksp,&nonzero);
69: if (nonzero) {
70: VecCopy(x,eis->b[pc->presolvedone-1]);
71: MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
72: }
74: /* save true b, other option is to swap pointers */
75: VecCopy(b,eis->b[pc->presolvedone-1]);
77: /* modify b by (L + D/omega)^{-1} */
78: MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP),0.0,1,1,b);
79: return(0);
80: }
82: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
83: {
84: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
88: /* get back true b */
89: VecCopy(eis->b[pc->presolvedone],b);
91: /* modify x by (U + D/omega)^{-1} */
92: VecCopy(x,eis->b[pc->presolvedone]);
93: MatSOR(eis->A,eis->b[pc->presolvedone],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP),0.0,1,1,x);
94: if (!pc->presolvedone) pc->mat = eis->A;
95: return(0);
96: }
98: static PetscErrorCode PCReset_Eisenstat(PC pc)
99: {
100: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
104: VecDestroy(&eis->b[0]);
105: VecDestroy(&eis->b[1]);
106: MatDestroy(&eis->shell);
107: VecDestroy(&eis->diag);
108: return(0);
109: }
111: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
112: {
116: PCReset_Eisenstat(pc);
117: PetscFree(pc->data);
118: return(0);
119: }
121: static PetscErrorCode PCSetFromOptions_Eisenstat(PetscOptionItems *PetscOptionsObject,PC pc)
122: {
123: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
125: PetscBool set,flg;
128: PetscOptionsHead(PetscOptionsObject,"Eisenstat SSOR options");
129: PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,NULL);
130: PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatSetNoDiagonalScaling",eis->usediag ? PETSC_FALSE : PETSC_TRUE,&flg,&set);
131: if (set) {
132: PCEisenstatSetNoDiagonalScaling(pc,flg);
133: }
134: PetscOptionsTail();
135: return(0);
136: }
138: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
139: {
140: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
142: PetscBool iascii;
145: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
146: if (iascii) {
147: PetscViewerASCIIPrintf(viewer," omega = %g\n",(double)eis->omega);
148: if (eis->usediag) {
149: PetscViewerASCIIPrintf(viewer," Using diagonal scaling (default)\n");
150: } else {
151: PetscViewerASCIIPrintf(viewer," Not using diagonal scaling\n");
152: }
153: }
154: return(0);
155: }
157: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
158: {
160: PetscInt M,N,m,n;
161: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
164: if (!pc->setupcalled) {
165: MatGetSize(pc->mat,&M,&N);
166: MatGetLocalSize(pc->mat,&m,&n);
167: MatCreate(PetscObjectComm((PetscObject)pc),&eis->shell);
168: MatSetSizes(eis->shell,m,n,M,N);
169: MatSetType(eis->shell,MATSHELL);
170: MatSetUp(eis->shell);
171: MatShellSetContext(eis->shell,pc);
172: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->shell);
173: MatShellSetOperation(eis->shell,MATOP_MULT,(void (*)(void))PCMult_Eisenstat);
174: }
175: if (!eis->usediag) return(0);
176: if (!pc->setupcalled) {
177: MatCreateVecs(pc->pmat,&eis->diag,NULL);
178: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->diag);
179: }
180: MatGetDiagonal(pc->pmat,eis->diag);
181: return(0);
182: }
184: /* --------------------------------------------------------------------*/
186: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
187: {
188: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
191: if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
192: eis->omega = omega;
193: return(0);
194: }
196: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc,PetscBool flg)
197: {
198: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
201: eis->usediag = flg;
202: return(0);
203: }
205: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc,PetscReal *omega)
206: {
207: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
210: *omega = eis->omega;
211: return(0);
212: }
214: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc,PetscBool *flg)
215: {
216: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
219: *flg = eis->usediag;
220: return(0);
221: }
223: /*@
224: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
225: to use with Eisenstat's trick (where omega = 1.0 by default).
227: Logically Collective on PC
229: Input Parameters:
230: + pc - the preconditioner context
231: - omega - relaxation coefficient (0 < omega < 2)
233: Options Database Key:
234: . -pc_eisenstat_omega <omega> - Sets omega
236: Notes:
237: The Eisenstat trick implementation of SSOR requires about 50% of the
238: usual amount of floating point operations used for SSOR + Krylov method;
239: however, the preconditioned problem must be solved with both left
240: and right preconditioning.
242: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
243: which can be chosen with the database options
244: $ -pc_type sor -pc_sor_symmetric
246: Level: intermediate
248: .seealso: PCSORSetOmega()
249: @*/
250: PetscErrorCode PCEisenstatSetOmega(PC pc,PetscReal omega)
251: {
257: PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
258: return(0);
259: }
261: /*@
262: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner
263: not to do additional diagonal preconditioning. For matrices with a constant
264: along the diagonal, this may save a small amount of work.
266: Logically Collective on PC
268: Input Parameters:
269: + pc - the preconditioner context
270: - flg - PETSC_TRUE turns off diagonal scaling inside the algorithm
272: Options Database Key:
273: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
275: Level: intermediate
277: Note:
278: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
279: likley want to use this routine since it will save you some unneeded flops.
281: .seealso: PCEisenstatSetOmega()
282: @*/
283: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)
284: {
289: PetscTryMethod(pc,"PCEisenstatSetNoDiagonalScaling_C",(PC,PetscBool),(pc,flg));
290: return(0);
291: }
293: /*@
294: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
295: to use with Eisenstat's trick (where omega = 1.0 by default).
297: Logically Collective on PC
299: Input Parameter:
300: . pc - the preconditioner context
302: Output Parameter:
303: . omega - relaxation coefficient (0 < omega < 2)
305: Options Database Key:
306: . -pc_eisenstat_omega <omega> - Sets omega
308: Notes:
309: The Eisenstat trick implementation of SSOR requires about 50% of the
310: usual amount of floating point operations used for SSOR + Krylov method;
311: however, the preconditioned problem must be solved with both left
312: and right preconditioning.
314: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
315: which can be chosen with the database options
316: $ -pc_type sor -pc_sor_symmetric
318: Level: intermediate
320: .seealso: PCSORGetOmega(), PCEisenstatSetOmega()
321: @*/
322: PetscErrorCode PCEisenstatGetOmega(PC pc,PetscReal *omega)
323: {
328: PetscUseMethod(pc,"PCEisenstatGetOmega_C",(PC,PetscReal*),(pc,omega));
329: return(0);
330: }
332: /*@
333: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
334: not to do additional diagonal preconditioning. For matrices with a constant
335: along the diagonal, this may save a small amount of work.
337: Logically Collective on PC
339: Input Parameter:
340: . pc - the preconditioner context
342: Output Parameter:
343: . flg - PETSC_TRUE means there is no diagonal scaling applied
345: Options Database Key:
346: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
348: Level: intermediate
350: Note:
351: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
352: likley want to use this routine since it will save you some unneeded flops.
354: .seealso: PCEisenstatGetOmega()
355: @*/
356: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool *flg)
357: {
362: PetscUseMethod(pc,"PCEisenstatGetNoDiagonalScaling_C",(PC,PetscBool*),(pc,flg));
363: return(0);
364: }
366: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool* change)
367: {
369: *change = PETSC_TRUE;
370: return(0);
371: }
373: /* --------------------------------------------------------------------*/
375: /*MC
376: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
377: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
379: Options Database Keys:
380: + -pc_eisenstat_omega <omega> - Sets omega
381: - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
383: Level: beginner
385: Notes:
386: Only implemented for the SeqAIJ matrix format.
387: Not a true parallel SOR, in parallel this implementation corresponds to block
388: Jacobi with SOR on each block.
390: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
391: PCEisenstatSetNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
392: M*/
394: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
395: {
397: PC_Eisenstat *eis;
400: PetscNewLog(pc,&eis);
402: pc->ops->apply = PCApply_Eisenstat;
403: pc->ops->presolve = PCPreSolve_Eisenstat;
404: pc->ops->postsolve = PCPostSolve_Eisenstat;
405: pc->ops->applyrichardson = NULL;
406: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
407: pc->ops->destroy = PCDestroy_Eisenstat;
408: pc->ops->reset = PCReset_Eisenstat;
409: pc->ops->view = PCView_Eisenstat;
410: pc->ops->setup = PCSetUp_Eisenstat;
412: pc->data = eis;
413: eis->omega = 1.0;
414: eis->b[0] = NULL;
415: eis->b[1] = NULL;
416: eis->diag = NULL;
417: eis->usediag = PETSC_TRUE;
419: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetOmega_C",PCEisenstatSetOmega_Eisenstat);
420: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetNoDiagonalScaling_C",PCEisenstatSetNoDiagonalScaling_Eisenstat);
421: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetOmega_C",PCEisenstatGetOmega_Eisenstat);
422: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetNoDiagonalScaling_C",PCEisenstatGetNoDiagonalScaling_Eisenstat);
423: PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_Eisenstat);
424: return(0);
425: }