Actual source code: eisen.c
1: /*$Id: eisen.c,v 1.113 2001/08/06 21:16:28 bsmith Exp $*/
3: /*
4: Defines a Eisenstat trick SSOR preconditioner. This uses about
5: %50 of the usual amount of floating point ops used for SSOR + Krylov
6: method. But it requires actually solving the preconditioned problem
7: with both left and right preconditioning.
8: */
9: #include src/sles/pc/pcimpl.h
11: typedef struct {
12: Mat shell,A;
13: Vec b,diag; /* temporary storage for true right hand side */
14: PetscReal omega;
15: PetscTruth usediag; /* indicates preconditioner should include diagonal scaling*/
16: } PC_Eisenstat;
19: #undef __FUNCT__
21: static int PCMult_Eisenstat(Mat mat,Vec b,Vec x)
22: {
23: int ierr;
24: PC pc;
25: PC_Eisenstat *eis;
28: MatShellGetContext(mat,(void **)&pc);
29: eis = (PC_Eisenstat*)pc->data;
30: MatRelax(eis->A,b,eis->omega,SOR_EISENSTAT,0.0,1,1,x);
31: return(0);
32: }
34: #undef __FUNCT__
36: static int PCApply_Eisenstat(PC pc,Vec x,Vec y)
37: {
38: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
39: int ierr;
42: if (eis->usediag) {VecPointwiseMult(x,eis->diag,y);}
43: else {VecCopy(x,y);}
44: return(0);
45: }
47: #undef __FUNCT__
49: static int PCPre_Eisenstat(PC pc,KSP ksp,Vec x,Vec b)
50: {
51: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
52: PetscTruth nonzero;
53: int ierr;
56: if (pc->mat != pc->pmat) SETERRQ(PETSC_ERR_SUP,"Cannot have different mat and pmat");
57:
58: /* swap shell matrix and true matrix */
59: eis->A = pc->mat;
60: pc->mat = eis->shell;
62: if (!eis->b) {
63: VecDuplicate(b,&eis->b);
64: PetscLogObjectParent(pc,eis->b);
65: }
66:
67: /* save true b, other option is to swap pointers */
68: VecCopy(b,eis->b);
70: /* if nonzero initial guess, modify x */
71: KSPGetInitialGuessNonzero(ksp,&nonzero);
72: if (nonzero) {
73: MatRelax(eis->A,x,eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
74: }
76: /* modify b by (L + D)^{-1} */
77: MatRelax(eis->A,b,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS |
78: SOR_FORWARD_SWEEP),0.0,1,1,b);
79: return(0);
80: }
82: #undef __FUNCT__
84: static int PCPost_Eisenstat(PC pc,KSP ksp,Vec x,Vec b)
85: {
86: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
87: int ierr;
90: MatRelax(eis->A,x,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS |
91: SOR_BACKWARD_SWEEP),0.0,1,1,x);
92: pc->mat = eis->A;
93: /* get back true b */
94: VecCopy(eis->b,b);
95: return(0);
96: }
98: #undef __FUNCT__
100: static int PCDestroy_Eisenstat(PC pc)
101: {
102: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
103: int ierr;
106: if (eis->b) {VecDestroy(eis->b);}
107: if (eis->shell) {MatDestroy(eis->shell);}
108: if (eis->diag) {VecDestroy(eis->diag);}
109: PetscFree(eis);
110: return(0);
111: }
113: #undef __FUNCT__
115: static int PCSetFromOptions_Eisenstat(PC pc)
116: {
117: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
118: int ierr;
119: PetscTruth flg;
122: PetscOptionsHead("Eisenstat SSOR options");
123: PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,0);
124: PetscOptionsName("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatNoDiagonalScaling",&flg);
125: if (flg) {
126: PCEisenstatNoDiagonalScaling(pc);
127: }
128: PetscOptionsTail();
129: return(0);
130: }
132: #undef __FUNCT__
134: static int PCView_Eisenstat(PC pc,PetscViewer viewer)
135: {
136: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
137: int ierr;
138: PetscTruth isascii;
141: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
142: if (isascii) {
143: PetscViewerASCIIPrintf(viewer,"Eisenstat: omega = %gn",eis->omega);
144: if (eis->usediag) {
145: PetscViewerASCIIPrintf(viewer,"Eisenstat: Using diagonal scaling (default)n");
146: } else {
147: PetscViewerASCIIPrintf(viewer,"Eisenstat: Not using diagonal scalingn");
148: }
149: } else {
150: SETERRQ1(1,"Viewer type %s not supported for Eisenstat PC",((PetscObject)viewer)->type_name);
151: }
152: return(0);
153: }
155: #undef __FUNCT__
157: static int PCSetUp_Eisenstat(PC pc)
158: {
159: int ierr,M,N,m,n;
160: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
163: if (!pc->setupcalled) {
164: MatGetSize(pc->mat,&M,&N);
165: MatGetLocalSize(pc->mat,&m,&n);
166: MatCreateShell(pc->comm,m,N,M,N,(void*)pc,&eis->shell);
167: PetscLogObjectParent(pc,eis->shell);
168: MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);
169: }
170: if (!eis->usediag) return(0);
171: if (!pc->setupcalled) {
172: VecDuplicate(pc->vec,&eis->diag);
173: PetscLogObjectParent(pc,eis->diag);
174: }
175: MatGetDiagonal(pc->pmat,eis->diag);
176: return(0);
177: }
179: /* --------------------------------------------------------------------*/
181: EXTERN_C_BEGIN
182: #undef __FUNCT__
184: int PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
185: {
186: PC_Eisenstat *eis;
189: if (omega >= 2.0 || omega <= 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
190: eis = (PC_Eisenstat*)pc->data;
191: eis->omega = omega;
192: return(0);
193: }
194: EXTERN_C_END
196: EXTERN_C_BEGIN
197: #undef __FUNCT__
199: int PCEisenstatNoDiagonalScaling_Eisenstat(PC pc)
200: {
201: PC_Eisenstat *eis;
204: eis = (PC_Eisenstat*)pc->data;
205: eis->usediag = PETSC_FALSE;
206: return(0);
207: }
208: EXTERN_C_END
210: #undef __FUNCT__
212: /*@
213: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
214: to use with Eisenstat's trick (where omega = 1.0 by default).
216: Collective on PC
218: Input Parameters:
219: + pc - the preconditioner context
220: - omega - relaxation coefficient (0 < omega < 2)
222: Options Database Key:
223: . -pc_eisenstat_omega <omega> - Sets omega
225: Notes:
226: The Eisenstat trick implementation of SSOR requires about 50% of the
227: usual amount of floating point operations used for SSOR + Krylov method;
228: however, the preconditioned problem must be solved with both left
229: and right preconditioning.
231: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
232: which can be chosen with the database options
233: $ -pc_type sor -pc_sor_symmetric
235: Level: intermediate
237: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega
239: .seealso: PCSORSetOmega()
240: @*/
241: int PCEisenstatSetOmega(PC pc,PetscReal omega)
242: {
243: int ierr,(*f)(PC,PetscReal);
247: PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatSetOmega_C",(void (**)(void))&f);
248: if (f) {
249: (*f)(pc,omega);
250: }
251: return(0);
252: }
254: #undef __FUNCT__
256: /*@
257: PCEisenstatNoDiagonalScaling - Causes the Eisenstat preconditioner
258: not to do additional diagonal preconditioning. For matrices with a constant
259: along the diagonal, this may save a small amount of work.
261: Collective on PC
263: Input Parameter:
264: . pc - the preconditioner context
266: Options Database Key:
267: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()
269: Level: intermediate
271: Note:
272: If you use the SLESSetDiagonalScaling() or -sles_diagonal_scale option then you will
273: likley want to use this routine since it will save you some unneeded flops.
275: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR
277: .seealso: PCEisenstatSetOmega()
278: @*/
279: int PCEisenstatNoDiagonalScaling(PC pc)
280: {
281: int ierr,(*f)(PC);
285: PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",(void (**)(void))&f);
286: if (f) {
287: (*f)(pc);
288: }
289: return(0);
290: }
292: /* --------------------------------------------------------------------*/
294: /*MC
295: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
296: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
298: Options Database Keys:
299: + -pc_eisenstat_omega <omega> - Sets omega
300: - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()
302: Level: beginner
304: Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick
306: Notes: Only implemented for the SeqAIJ matrix format.
307: Not a true parallel SOR, in parallel this implementation corresponds to block
308: Jacobi with SOR on each block.
310: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
311: PCEisenstatNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
312: M*/
314: EXTERN_C_BEGIN
315: #undef __FUNCT__
317: int PCCreate_Eisenstat(PC pc)
318: {
319: int ierr;
320: PC_Eisenstat *eis;
323: PetscNew(PC_Eisenstat,&eis);
324: PetscLogObjectMemory(pc,sizeof(PC_Eisenstat));
326: pc->ops->apply = PCApply_Eisenstat;
327: pc->ops->presolve = PCPre_Eisenstat;
328: pc->ops->postsolve = PCPost_Eisenstat;
329: pc->ops->applyrichardson = 0;
330: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
331: pc->ops->destroy = PCDestroy_Eisenstat;
332: pc->ops->view = PCView_Eisenstat;
333: pc->ops->setup = PCSetUp_Eisenstat;
335: pc->data = (void*)eis;
336: eis->omega = 1.0;
337: eis->b = 0;
338: eis->diag = 0;
339: eis->usediag = PETSC_TRUE;
341: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatSetOmega_C","PCEisenstatSetOmega_Eisenstat",
342: PCEisenstatSetOmega_Eisenstat);
343: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",
344: "PCEisenstatNoDiagonalScaling_Eisenstat",
345: PCEisenstatNoDiagonalScaling_Eisenstat);
346: return(0);
347: }
348: EXTERN_C_END