Actual source code: eisen.c
petsc-3.6.1 2015-08-06
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(PetscObjectComm((PetscObject)pc),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((PetscObject)pc,(PetscObject)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) pc->mat = eis->A;
104: return(0);
105: }
109: static PetscErrorCode PCReset_Eisenstat(PC pc)
110: {
111: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
115: VecDestroy(&eis->b[0]);
116: VecDestroy(&eis->b[1]);
117: MatDestroy(&eis->shell);
118: VecDestroy(&eis->diag);
119: return(0);
120: }
124: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
125: {
129: PCReset_Eisenstat(pc);
130: PetscFree(pc->data);
131: return(0);
132: }
136: static PetscErrorCode PCSetFromOptions_Eisenstat(PetscOptions *PetscOptionsObject,PC pc)
137: {
138: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
140: PetscBool set,flg;
143: PetscOptionsHead(PetscOptionsObject,"Eisenstat SSOR options");
144: PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,NULL);
145: PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatSetNoDiagonalScaling",eis->usediag ? PETSC_FALSE : PETSC_TRUE,&flg,&set);
146: if (set) {
147: PCEisenstatSetNoDiagonalScaling(pc,flg);
148: }
149: PetscOptionsTail();
150: return(0);
151: }
155: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
156: {
157: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
159: PetscBool iascii;
162: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
163: if (iascii) {
164: PetscViewerASCIIPrintf(viewer,"Eisenstat: omega = %g\n",(double)eis->omega);
165: if (eis->usediag) {
166: PetscViewerASCIIPrintf(viewer,"Eisenstat: Using diagonal scaling (default)\n");
167: } else {
168: PetscViewerASCIIPrintf(viewer,"Eisenstat: Not using diagonal scaling\n");
169: }
170: }
171: return(0);
172: }
176: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
177: {
179: PetscInt M,N,m,n;
180: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
183: if (!pc->setupcalled) {
184: MatGetSize(pc->mat,&M,&N);
185: MatGetLocalSize(pc->mat,&m,&n);
186: MatCreate(PetscObjectComm((PetscObject)pc),&eis->shell);
187: MatSetSizes(eis->shell,m,n,M,N);
188: MatSetType(eis->shell,MATSHELL);
189: MatSetUp(eis->shell);
190: MatShellSetContext(eis->shell,(void*)pc);
191: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->shell);
192: MatShellSetOperation(eis->shell,MATOP_MULT,(void (*)(void))PCMult_Eisenstat);
193: }
194: if (!eis->usediag) return(0);
195: if (!pc->setupcalled) {
196: MatCreateVecs(pc->pmat,&eis->diag,0);
197: PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->diag);
198: }
199: MatGetDiagonal(pc->pmat,eis->diag);
200: return(0);
201: }
203: /* --------------------------------------------------------------------*/
207: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
208: {
209: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
212: if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
213: eis->omega = omega;
214: return(0);
215: }
219: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc,PetscBool flg)
220: {
221: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
224: eis->usediag = flg;
225: return(0);
226: }
230: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc,PetscReal *omega)
231: {
232: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
235: *omega = eis->omega;
236: return(0);
237: }
241: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc,PetscBool *flg)
242: {
243: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
246: *flg = eis->usediag;
247: return(0);
248: }
252: /*@
253: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
254: to use with Eisenstat's trick (where omega = 1.0 by default).
256: Logically Collective on PC
258: Input Parameters:
259: + pc - the preconditioner context
260: - omega - relaxation coefficient (0 < omega < 2)
262: Options Database Key:
263: . -pc_eisenstat_omega <omega> - Sets omega
265: Notes:
266: The Eisenstat trick implementation of SSOR requires about 50% of the
267: usual amount of floating point operations used for SSOR + Krylov method;
268: however, the preconditioned problem must be solved with both left
269: and right preconditioning.
271: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
272: which can be chosen with the database options
273: $ -pc_type sor -pc_sor_symmetric
275: Level: intermediate
277: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega
279: .seealso: PCSORSetOmega()
280: @*/
281: PetscErrorCode PCEisenstatSetOmega(PC pc,PetscReal omega)
282: {
288: PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
289: return(0);
290: }
294: /*@
295: PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner
296: not to do additional diagonal preconditioning. For matrices with a constant
297: along the diagonal, this may save a small amount of work.
299: Logically Collective on PC
301: Input Parameters:
302: + pc - the preconditioner context
303: - flg - PETSC_TRUE turns off diagonal scaling inside the algorithm
305: Options Database Key:
306: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
308: Level: intermediate
310: Note:
311: If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
312: likley want to use this routine since it will save you some unneeded flops.
314: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR
316: .seealso: PCEisenstatSetOmega()
317: @*/
318: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)
319: {
324: PetscTryMethod(pc,"PCEisenstatSetNoDiagonalScaling_C",(PC,PetscBool),(pc,flg));
325: return(0);
326: }
330: /*@
331: PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
332: to use with Eisenstat's trick (where omega = 1.0 by default).
334: Logically Collective on PC
336: Input Parameter:
337: . pc - the preconditioner context
339: Output Parameter:
340: . omega - relaxation coefficient (0 < omega < 2)
342: Options Database Key:
343: . -pc_eisenstat_omega <omega> - Sets omega
345: Notes:
346: The Eisenstat trick implementation of SSOR requires about 50% of the
347: usual amount of floating point operations used for SSOR + Krylov method;
348: however, the preconditioned problem must be solved with both left
349: and right preconditioning.
351: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
352: which can be chosen with the database options
353: $ -pc_type sor -pc_sor_symmetric
355: Level: intermediate
357: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega
359: .seealso: PCSORGetOmega(), PCEisenstatSetOmega()
360: @*/
361: PetscErrorCode PCEisenstatGetOmega(PC pc,PetscReal *omega)
362: {
367: PetscUseMethod(pc,"PCEisenstatGetOmega_C",(PC,PetscReal*),(pc,omega));
368: return(0);
369: }
373: /*@
374: PCEisenstatGetNoDiagonalScaling - Causes the Eisenstat preconditioner
375: not to do additional diagonal preconditioning. For matrices with a constant
376: along the diagonal, this may save a small amount of work.
378: Logically Collective on PC
380: Input Parameter:
381: . pc - the preconditioner context
383: Output Parameter:
384: . flg - PETSC_TRUE means there is no diagonal scaling applied
386: Options Database Key:
387: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
389: Level: intermediate
391: Note:
392: If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
393: likley want to use this routine since it will save you some unneeded flops.
395: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR
397: .seealso: PCEisenstatGetOmega()
398: @*/
399: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool *flg)
400: {
405: PetscTryMethod(pc,"PCEisenstatGetNoDiagonalScaling_C",(PC,PetscBool*),(pc,flg));
406: return(0);
407: }
409: /* --------------------------------------------------------------------*/
411: /*MC
412: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
413: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
415: Options Database Keys:
416: + -pc_eisenstat_omega <omega> - Sets omega
417: - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()
419: Level: beginner
421: Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick
423: Notes: Only implemented for the SeqAIJ matrix format.
424: Not a true parallel SOR, in parallel this implementation corresponds to block
425: Jacobi with SOR on each block.
427: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
428: PCEisenstatSetNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
429: M*/
433: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
434: {
436: PC_Eisenstat *eis;
439: PetscNewLog(pc,&eis);
441: pc->ops->apply = PCApply_Eisenstat;
442: pc->ops->presolve = PCPreSolve_Eisenstat;
443: pc->ops->postsolve = PCPostSolve_Eisenstat;
444: pc->ops->applyrichardson = 0;
445: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
446: pc->ops->destroy = PCDestroy_Eisenstat;
447: pc->ops->reset = PCReset_Eisenstat;
448: pc->ops->view = PCView_Eisenstat;
449: pc->ops->setup = PCSetUp_Eisenstat;
451: pc->data = (void*)eis;
452: eis->omega = 1.0;
453: eis->b[0] = 0;
454: eis->b[1] = 0;
455: eis->diag = 0;
456: eis->usediag = PETSC_TRUE;
458: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetOmega_C",PCEisenstatSetOmega_Eisenstat);
459: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetNoDiagonalScaling_C",PCEisenstatSetNoDiagonalScaling_Eisenstat);
460: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetOmega_C",PCEisenstatGetOmega_Eisenstat);
461: PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetNoDiagonalScaling_C",PCEisenstatGetNoDiagonalScaling_Eisenstat);
462: return(0);
463: }