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;
18: static PetscErrorCode PCMult_Eisenstat(Mat mat,Vec b,Vec x)
19: {
21: PC pc;
22: PC_Eisenstat *eis;
25: MatShellGetContext(mat,(void**)&pc);
26: eis = (PC_Eisenstat*)pc->data;
27: MatSOR(eis->A,b,eis->omega,SOR_EISENSTAT,0.0,1,1,x);
28: return(0);
29: }
31: static PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y)
32: {
33: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
35: PetscBool hasop;
38: if (eis->usediag) {
39: MatHasOperation(pc->pmat,MATOP_MULT_DIAGONAL_BLOCK,&hasop);
40: if (hasop) {
41: MatMultDiagonalBlock(pc->pmat,x,y);
42: } else {
43: VecPointwiseMult(y,x,eis->diag);
44: }
45: } else {VecCopy(x,y);}
46: return(0);
47: }
49: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
50: {
51: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
52: PetscBool nonzero;
56: if (pc->presolvedone < 2) {
57: if (pc->mat != pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot have different mat and pmat");
58: /* swap shell matrix and true matrix */
59: eis->A = pc->mat;
60: pc->mat = eis->shell;
61: }
63: if (!eis->b[pc->presolvedone-1]) {
64: VecDuplicate(b,&eis->b[pc->presolvedone-1]);
65: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->b[pc->presolvedone-1]);
66: }
68: /* if nonzero initial guess, modify x */
69: KSPGetInitialGuessNonzero(ksp,&nonzero);
70: if (nonzero) {
71: VecCopy(x,eis->b[pc->presolvedone-1]);
72: MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
73: }
75: /* save true b, other option is to swap pointers */
76: VecCopy(b,eis->b[pc->presolvedone-1]);
78: /* modify b by (L + D/omega)^{-1} */
79: MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP),0.0,1,1,b);
80: return(0);
81: }
83: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
84: {
85: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
89: /* get back true b */
90: VecCopy(eis->b[pc->presolvedone],b);
92: /* modify x by (U + D/omega)^{-1} */
93: VecCopy(x,eis->b[pc->presolvedone]);
94: MatSOR(eis->A,eis->b[pc->presolvedone],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP),0.0,1,1,x);
95: if (!pc->presolvedone) pc->mat = eis->A;
96: return(0);
97: }
99: static PetscErrorCode PCReset_Eisenstat(PC pc)
100: {
101: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
105: VecDestroy(&eis->b[0]);
106: VecDestroy(&eis->b[1]);
107: MatDestroy(&eis->shell);
108: VecDestroy(&eis->diag);
109: return(0);
110: }
112: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
113: {
117: PCReset_Eisenstat(pc);
118: PetscFree(pc->data);
119: return(0);
120: }
122: static PetscErrorCode PCSetFromOptions_Eisenstat(PetscOptionItems *PetscOptionsObject,PC pc)
123: {
124: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
126: PetscBool set,flg;
129: PetscOptionsHead(PetscOptionsObject,"Eisenstat SSOR options");
130: PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,NULL);
131: PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatSetNoDiagonalScaling",eis->usediag ? PETSC_FALSE : PETSC_TRUE,&flg,&set);
132: if (set) {
133: PCEisenstatSetNoDiagonalScaling(pc,flg);
134: }
135: PetscOptionsTail();
136: return(0);
137: }
139: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
140: {
141: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
143: PetscBool iascii;
146: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
147: if (iascii) {
148: PetscViewerASCIIPrintf(viewer," omega = %g\n",(double)eis->omega);
149: if (eis->usediag) {
150: PetscViewerASCIIPrintf(viewer," Using diagonal scaling (default)\n");
151: } else {
152: PetscViewerASCIIPrintf(viewer," Not using diagonal scaling\n");
153: }
154: }
155: return(0);
156: }
158: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
159: {
161: PetscInt M,N,m,n;
162: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
165: if (!pc->setupcalled) {
166: MatGetSize(pc->mat,&M,&N);
167: MatGetLocalSize(pc->mat,&m,&n);
168: MatCreate(PetscObjectComm((PetscObject)pc),&eis->shell);
169: MatSetSizes(eis->shell,m,n,M,N);
170: MatSetType(eis->shell,MATSHELL);
171: MatSetUp(eis->shell);
172: MatShellSetContext(eis->shell,(void*)pc);
173: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->shell);
174: MatShellSetOperation(eis->shell,MATOP_MULT,(void (*)(void))PCMult_Eisenstat);
175: }
176: if (!eis->usediag) return(0);
177: if (!pc->setupcalled) {
178: MatCreateVecs(pc->pmat,&eis->diag,NULL);
179: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->diag);
180: }
181: MatGetDiagonal(pc->pmat,eis->diag);
182: return(0);
183: }
185: /* --------------------------------------------------------------------*/
187: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
188: {
189: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
192: if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
193: eis->omega = omega;
194: return(0);
195: }
197: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc,PetscBool flg)
198: {
199: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
202: eis->usediag = flg;
203: return(0);
204: }
206: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc,PetscReal *omega)
207: {
208: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
211: *omega = eis->omega;
212: return(0);
213: }
215: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc,PetscBool *flg)
216: {
217: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
220: *flg = eis->usediag;
221: return(0);
222: }
224: /*@
225: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
226: to use with Eisenstat's trick (where omega = 1.0 by default).
228: Logically Collective on PC
230: Input Parameters:
231: + pc - the preconditioner context
232: - omega - relaxation coefficient (0 < omega < 2)
234: Options Database Key:
235: . -pc_eisenstat_omega <omega> - Sets omega
237: Notes:
238: The Eisenstat trick implementation of SSOR requires about 50% of the
239: usual amount of floating point operations used for SSOR + Krylov method;
240: however, the preconditioned problem must be solved with both left
241: and right preconditioning.
243: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
244: which can be chosen with the database options
245: $ -pc_type sor -pc_sor_symmetric
247: Level: intermediate
249: .seealso: PCSORSetOmega()
250: @*/
251: PetscErrorCode PCEisenstatSetOmega(PC pc,PetscReal omega)
252: {
258: PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
259: return(0);
260: }
262: /*@
263: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner
264: not to do additional diagonal preconditioning. For matrices with a constant
265: along the diagonal, this may save a small amount of work.
267: Logically Collective on PC
269: Input Parameters:
270: + pc - the preconditioner context
271: - flg - PETSC_TRUE turns off diagonal scaling inside the algorithm
273: Options Database Key:
274: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
276: Level: intermediate
278: Note:
279: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
280: likley want to use this routine since it will save you some unneeded flops.
282: .seealso: PCEisenstatSetOmega()
283: @*/
284: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)
285: {
290: PetscTryMethod(pc,"PCEisenstatSetNoDiagonalScaling_C",(PC,PetscBool),(pc,flg));
291: return(0);
292: }
294: /*@
295: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
296: to use with Eisenstat's trick (where omega = 1.0 by default).
298: Logically Collective on PC
300: Input Parameter:
301: . pc - the preconditioner context
303: Output Parameter:
304: . omega - relaxation coefficient (0 < omega < 2)
306: Options Database Key:
307: . -pc_eisenstat_omega <omega> - Sets omega
309: Notes:
310: The Eisenstat trick implementation of SSOR requires about 50% of the
311: usual amount of floating point operations used for SSOR + Krylov method;
312: however, the preconditioned problem must be solved with both left
313: and right preconditioning.
315: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
316: which can be chosen with the database options
317: $ -pc_type sor -pc_sor_symmetric
319: Level: intermediate
321: .seealso: PCSORGetOmega(), PCEisenstatSetOmega()
322: @*/
323: PetscErrorCode PCEisenstatGetOmega(PC pc,PetscReal *omega)
324: {
329: PetscUseMethod(pc,"PCEisenstatGetOmega_C",(PC,PetscReal*),(pc,omega));
330: return(0);
331: }
333: /*@
334: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
335: not to do additional diagonal preconditioning. For matrices with a constant
336: along the diagonal, this may save a small amount of work.
338: Logically Collective on PC
340: Input Parameter:
341: . pc - the preconditioner context
343: Output Parameter:
344: . flg - PETSC_TRUE means there is no diagonal scaling applied
346: Options Database Key:
347: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
349: Level: intermediate
351: Note:
352: If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
353: likley want to use this routine since it will save you some unneeded flops.
355: .seealso: PCEisenstatGetOmega()
356: @*/
357: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool *flg)
358: {
363: PetscUseMethod(pc,"PCEisenstatGetNoDiagonalScaling_C",(PC,PetscBool*),(pc,flg));
364: return(0);
365: }
367: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool* change)
368: {
370: *change = PETSC_TRUE;
371: return(0);
372: }
374: /* --------------------------------------------------------------------*/
376: /*MC
377: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
378: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
380: Options Database Keys:
381: + -pc_eisenstat_omega <omega> - Sets omega
382: - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
384: Level: beginner
386: Notes:
387: Only implemented for the SeqAIJ matrix format.
388: Not a true parallel SOR, in parallel this implementation corresponds to block
389: Jacobi with SOR on each block.
391: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
392: PCEisenstatSetNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
393: M*/
395: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
396: {
398: PC_Eisenstat *eis;
401: PetscNewLog(pc,&eis);
403: pc->ops->apply = PCApply_Eisenstat;
404: pc->ops->presolve = PCPreSolve_Eisenstat;
405: pc->ops->postsolve = PCPostSolve_Eisenstat;
406: pc->ops->applyrichardson = NULL;
407: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
408: pc->ops->destroy = PCDestroy_Eisenstat;
409: pc->ops->reset = PCReset_Eisenstat;
410: pc->ops->view = PCView_Eisenstat;
411: pc->ops->setup = PCSetUp_Eisenstat;
413: pc->data = (void*)eis;
414: eis->omega = 1.0;
415: eis->b[0] = NULL;
416: eis->b[1] = NULL;
417: eis->diag = NULL;
418: eis->usediag = PETSC_TRUE;
420: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetOmega_C",PCEisenstatSetOmega_Eisenstat);
421: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetNoDiagonalScaling_C",PCEisenstatSetNoDiagonalScaling_Eisenstat);
422: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetOmega_C",PCEisenstatGetOmega_Eisenstat);
423: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetNoDiagonalScaling_C",PCEisenstatGetNoDiagonalScaling_Eisenstat);
424: PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_Eisenstat);
425: return(0);
426: }