Actual source code: eisen.c
petsc-3.9.4 2018-09-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>
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,0);
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: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega
251: .seealso: PCSORSetOmega()
252: @*/
253: PetscErrorCode PCEisenstatSetOmega(PC pc,PetscReal omega)
254: {
260: PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
261: return(0);
262: }
264: /*@
265: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner
266: not to do additional diagonal preconditioning. For matrices with a constant
267: along the diagonal, this may save a small amount of work.
269: Logically Collective on PC
271: Input Parameters:
272: + pc - the preconditioner context
273: - flg - PETSC_TRUE turns off diagonal scaling inside the algorithm
275: Options Database Key:
276: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
278: Level: intermediate
280: Note:
281: If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
282: likley want to use this routine since it will save you some unneeded flops.
284: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR
286: .seealso: PCEisenstatSetOmega()
287: @*/
288: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)
289: {
294: PetscTryMethod(pc,"PCEisenstatSetNoDiagonalScaling_C",(PC,PetscBool),(pc,flg));
295: return(0);
296: }
298: /*@
299: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
300: to use with Eisenstat's trick (where omega = 1.0 by default).
302: Logically Collective on PC
304: Input Parameter:
305: . pc - the preconditioner context
307: Output Parameter:
308: . omega - relaxation coefficient (0 < omega < 2)
310: Options Database Key:
311: . -pc_eisenstat_omega <omega> - Sets omega
313: Notes:
314: The Eisenstat trick implementation of SSOR requires about 50% of the
315: usual amount of floating point operations used for SSOR + Krylov method;
316: however, the preconditioned problem must be solved with both left
317: and right preconditioning.
319: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
320: which can be chosen with the database options
321: $ -pc_type sor -pc_sor_symmetric
323: Level: intermediate
325: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega
327: .seealso: PCSORGetOmega(), PCEisenstatSetOmega()
328: @*/
329: PetscErrorCode PCEisenstatGetOmega(PC pc,PetscReal *omega)
330: {
335: PetscUseMethod(pc,"PCEisenstatGetOmega_C",(PC,PetscReal*),(pc,omega));
336: return(0);
337: }
339: /*@
340: PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
341: not to do additional diagonal preconditioning. For matrices with a constant
342: along the diagonal, this may save a small amount of work.
344: Logically Collective on PC
346: Input Parameter:
347: . pc - the preconditioner context
349: Output Parameter:
350: . flg - PETSC_TRUE means there is no diagonal scaling applied
352: Options Database Key:
353: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
355: Level: intermediate
357: Note:
358: If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
359: likley want to use this routine since it will save you some unneeded flops.
361: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR
363: .seealso: PCEisenstatGetOmega()
364: @*/
365: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool *flg)
366: {
371: PetscUseMethod(pc,"PCEisenstatGetNoDiagonalScaling_C",(PC,PetscBool*),(pc,flg));
372: return(0);
373: }
375: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool* change)
376: {
378: *change = PETSC_TRUE;
379: return(0);
380: }
382: /* --------------------------------------------------------------------*/
384: /*MC
385: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
386: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
388: Options Database Keys:
389: + -pc_eisenstat_omega <omega> - Sets omega
390: - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
392: Level: beginner
394: Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick
396: Notes: Only implemented for the SeqAIJ matrix format.
397: Not a true parallel SOR, in parallel this implementation corresponds to block
398: Jacobi with SOR on each block.
400: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
401: PCEisenstatSetNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
402: M*/
404: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
405: {
407: PC_Eisenstat *eis;
410: PetscNewLog(pc,&eis);
412: pc->ops->apply = PCApply_Eisenstat;
413: pc->ops->presolve = PCPreSolve_Eisenstat;
414: pc->ops->postsolve = PCPostSolve_Eisenstat;
415: pc->ops->applyrichardson = 0;
416: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
417: pc->ops->destroy = PCDestroy_Eisenstat;
418: pc->ops->reset = PCReset_Eisenstat;
419: pc->ops->view = PCView_Eisenstat;
420: pc->ops->setup = PCSetUp_Eisenstat;
422: pc->data = (void*)eis;
423: eis->omega = 1.0;
424: eis->b[0] = 0;
425: eis->b[1] = 0;
426: eis->diag = 0;
427: eis->usediag = PETSC_TRUE;
429: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetOmega_C",PCEisenstatSetOmega_Eisenstat);
430: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetNoDiagonalScaling_C",PCEisenstatSetNoDiagonalScaling_Eisenstat);
431: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetOmega_C",PCEisenstatGetOmega_Eisenstat);
432: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetNoDiagonalScaling_C",PCEisenstatGetNoDiagonalScaling_Eisenstat);
433: PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_Eisenstat);
434: return(0);
435: }