Actual source code: eisen.c
petsc-3.3-p7 2013-05-11
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> /*I "petscpc.h" I*/
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;
20: static PetscErrorCode PCMult_Eisenstat(Mat mat,Vec b,Vec x)
21: {
23: PC pc;
24: PC_Eisenstat *eis;
27: MatShellGetContext(mat,(void **)&pc);
28: eis = (PC_Eisenstat*)pc->data;
29: MatSOR(eis->A,b,eis->omega,SOR_EISENSTAT,0.0,1,1,x);
30: return(0);
31: }
35: static PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y)
36: {
37: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
39: PetscBool hasop;
42: if (eis->usediag) {
43: MatHasOperation(pc->pmat,MATOP_MULT_DIAGONAL_BLOCK,&hasop);
44: if (hasop) {
45: MatMultDiagonalBlock(pc->pmat,x,y);
46: } else {
47: VecPointwiseMult(y,x,eis->diag);
48: }
49: } else {VecCopy(x,y);}
50: return(0);
51: }
55: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
56: {
57: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
58: PetscBool nonzero;
62: if (pc->presolvedone < 2) {
63: if (pc->mat != pc->pmat) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cannot have different mat and pmat");
64: /* swap shell matrix and true matrix */
65: eis->A = pc->mat;
66: pc->mat = eis->shell;
67: }
69: if (!eis->b[pc->presolvedone-1]) {
70: VecDuplicate(b,&eis->b[pc->presolvedone-1]);
71: PetscLogObjectParent(pc,eis->b[pc->presolvedone-1]);
72: }
74: /* if nonzero initial guess, modify x */
75: KSPGetInitialGuessNonzero(ksp,&nonzero);
76: if (nonzero) {
77: VecCopy(x,eis->b[pc->presolvedone-1]);
78: MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
79: }
81: /* save true b, other option is to swap pointers */
82: VecCopy(b,eis->b[pc->presolvedone-1]);
84: /* modify b by (L + D/omega)^{-1} */
85: MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP),0.0,1,1,b);
86: return(0);
87: }
91: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
92: {
93: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
97: /* get back true b */
98: VecCopy(eis->b[pc->presolvedone],b);
100: /* modify x by (U + D/omega)^{-1} */
101: VecCopy(x,eis->b[pc->presolvedone]);
102: MatSOR(eis->A,eis->b[pc->presolvedone],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP),0.0,1,1,x);
103: if (!pc->presolvedone) {
104: pc->mat = eis->A;
105: }
106: return(0);
107: }
111: static PetscErrorCode PCReset_Eisenstat(PC pc)
112: {
113: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
117: VecDestroy(&eis->b[0]);
118: VecDestroy(&eis->b[1]);
119: MatDestroy(&eis->shell);
120: VecDestroy(&eis->diag);
121: return(0);
122: }
126: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
127: {
131: PCReset_Eisenstat(pc);
132: PetscFree(pc->data);
133: return(0);
134: }
138: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc)
139: {
140: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
142: PetscBool flg = PETSC_FALSE;
145: PetscOptionsHead("Eisenstat SSOR options");
146: PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,PETSC_NULL);
147: PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatNoDiagonalScaling",flg,&flg,PETSC_NULL);
148: if (flg) {
149: PCEisenstatNoDiagonalScaling(pc);
150: }
151: PetscOptionsTail();
152: return(0);
153: }
157: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
158: {
159: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
161: PetscBool iascii;
164: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
165: if (iascii) {
166: PetscViewerASCIIPrintf(viewer,"Eisenstat: omega = %G\n",eis->omega);
167: if (eis->usediag) {
168: PetscViewerASCIIPrintf(viewer,"Eisenstat: Using diagonal scaling (default)\n");
169: } else {
170: PetscViewerASCIIPrintf(viewer,"Eisenstat: Not using diagonal scaling\n");
171: }
172: } else {
173: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for Eisenstat PC",((PetscObject)viewer)->type_name);
174: }
175: return(0);
176: }
180: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
181: {
183: PetscInt M,N,m,n;
184: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
187: if (!pc->setupcalled) {
188: MatGetSize(pc->mat,&M,&N);
189: MatGetLocalSize(pc->mat,&m,&n);
190: MatCreate(((PetscObject)pc)->comm,&eis->shell);
191: MatSetSizes(eis->shell,m,n,M,N);
192: MatSetType(eis->shell,MATSHELL);
193: MatSetUp(eis->shell);
194: MatShellSetContext(eis->shell,(void*)pc);
195: PetscLogObjectParent(pc,eis->shell);
196: MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);
197: }
198: if (!eis->usediag) return(0);
199: if (!pc->setupcalled) {
200: MatGetVecs(pc->pmat,&eis->diag,0);
201: PetscLogObjectParent(pc,eis->diag);
202: }
203: MatGetDiagonal(pc->pmat,eis->diag);
204: return(0);
205: }
207: /* --------------------------------------------------------------------*/
209: EXTERN_C_BEGIN
212: PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
213: {
214: PC_Eisenstat *eis;
217: if (omega >= 2.0 || omega <= 0.0) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
218: eis = (PC_Eisenstat*)pc->data;
219: eis->omega = omega;
220: return(0);
221: }
222: EXTERN_C_END
224: EXTERN_C_BEGIN
227: PetscErrorCode PCEisenstatNoDiagonalScaling_Eisenstat(PC pc)
228: {
229: PC_Eisenstat *eis;
232: eis = (PC_Eisenstat*)pc->data;
233: eis->usediag = PETSC_FALSE;
234: return(0);
235: }
236: EXTERN_C_END
240: /*@
241: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
242: to use with Eisenstat's trick (where omega = 1.0 by default).
244: Logically Collective on PC
246: Input Parameters:
247: + pc - the preconditioner context
248: - omega - relaxation coefficient (0 < omega < 2)
250: Options Database Key:
251: . -pc_eisenstat_omega <omega> - Sets omega
253: Notes:
254: The Eisenstat trick implementation of SSOR requires about 50% of the
255: usual amount of floating point operations used for SSOR + Krylov method;
256: however, the preconditioned problem must be solved with both left
257: and right preconditioning.
259: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
260: which can be chosen with the database options
261: $ -pc_type sor -pc_sor_symmetric
263: Level: intermediate
265: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega
267: .seealso: PCSORSetOmega()
268: @*/
269: PetscErrorCode PCEisenstatSetOmega(PC pc,PetscReal omega)
270: {
276: PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
277: return(0);
278: }
282: /*@
283: PCEisenstatNoDiagonalScaling - Causes the Eisenstat preconditioner
284: not to do additional diagonal preconditioning. For matrices with a constant
285: along the diagonal, this may save a small amount of work.
287: Logically Collective on PC
289: Input Parameter:
290: . pc - the preconditioner context
292: Options Database Key:
293: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()
295: Level: intermediate
297: Note:
298: If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
299: likley want to use this routine since it will save you some unneeded flops.
301: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR
303: .seealso: PCEisenstatSetOmega()
304: @*/
305: PetscErrorCode PCEisenstatNoDiagonalScaling(PC pc)
306: {
311: PetscTryMethod(pc,"PCEisenstatNoDiagonalScaling_C",(PC),(pc));
312: return(0);
313: }
315: /* --------------------------------------------------------------------*/
317: /*MC
318: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
319: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
321: Options Database Keys:
322: + -pc_eisenstat_omega <omega> - Sets omega
323: - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()
325: Level: beginner
327: Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick
329: Notes: Only implemented for the SeqAIJ matrix format.
330: Not a true parallel SOR, in parallel this implementation corresponds to block
331: Jacobi with SOR on each block.
333: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
334: PCEisenstatNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
335: M*/
337: EXTERN_C_BEGIN
340: PetscErrorCode PCCreate_Eisenstat(PC pc)
341: {
343: PC_Eisenstat *eis;
346: PetscNewLog(pc,PC_Eisenstat,&eis);
348: pc->ops->apply = PCApply_Eisenstat;
349: pc->ops->presolve = PCPreSolve_Eisenstat;
350: pc->ops->postsolve = PCPostSolve_Eisenstat;
351: pc->ops->applyrichardson = 0;
352: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
353: pc->ops->destroy = PCDestroy_Eisenstat;
354: pc->ops->reset = PCReset_Eisenstat;
355: pc->ops->view = PCView_Eisenstat;
356: pc->ops->setup = PCSetUp_Eisenstat;
358: pc->data = (void*)eis;
359: eis->omega = 1.0;
360: eis->b[0] = 0;
361: eis->b[1] = 0;
362: eis->diag = 0;
363: eis->usediag = PETSC_TRUE;
365: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatSetOmega_C","PCEisenstatSetOmega_Eisenstat",
366: PCEisenstatSetOmega_Eisenstat);
367: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",
368: "PCEisenstatNoDiagonalScaling_Eisenstat",
369: PCEisenstatNoDiagonalScaling_Eisenstat);
370: return(0);
371: }
372: EXTERN_C_END