Actual source code: eisen.c
1: #define PETSCKSP_DLL
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 private/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;
21: static PetscErrorCode PCMult_Eisenstat(Mat mat,Vec b,Vec x)
22: {
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: }
36: static PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y)
37: {
38: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
42: if (eis->usediag) {VecPointwiseMult(y,x,eis->diag);}
43: else {VecCopy(x,y);}
44: return(0);
45: }
49: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
50: {
51: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
52: PetscTruth nonzero;
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 | SOR_FORWARD_SWEEP),0.0,1,1,b);
78: return(0);
79: }
83: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
84: {
85: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
89: MatRelax(eis->A,x,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_BACKWARD_SWEEP),0.0,1,1,x);
90: pc->mat = eis->A;
91: /* get back true b */
92: VecCopy(eis->b,b);
93: return(0);
94: }
98: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
99: {
100: PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
104: if (eis->b) {VecDestroy(eis->b);}
105: if (eis->shell) {MatDestroy(eis->shell);}
106: if (eis->diag) {VecDestroy(eis->diag);}
107: PetscFree(eis);
108: return(0);
109: }
113: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc)
114: {
115: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
117: PetscTruth flg;
120: PetscOptionsHead("Eisenstat SSOR options");
121: PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,0);
122: PetscOptionsName("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatNoDiagonalScaling",&flg);
123: if (flg) {
124: PCEisenstatNoDiagonalScaling(pc);
125: }
126: PetscOptionsTail();
127: return(0);
128: }
132: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
133: {
134: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
136: PetscTruth iascii;
139: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
140: if (iascii) {
141: PetscViewerASCIIPrintf(viewer,"Eisenstat: omega = %G\n",eis->omega);
142: if (eis->usediag) {
143: PetscViewerASCIIPrintf(viewer,"Eisenstat: Using diagonal scaling (default)\n");
144: } else {
145: PetscViewerASCIIPrintf(viewer,"Eisenstat: Not using diagonal scaling\n");
146: }
147: } else {
148: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for Eisenstat PC",((PetscObject)viewer)->type_name);
149: }
150: return(0);
151: }
155: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
156: {
158: PetscInt M,N,m,n;
159: PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
162: if (!pc->setupcalled) {
163: MatGetSize(pc->mat,&M,&N);
164: MatGetLocalSize(pc->mat,&m,&n);
165: MatCreate(pc->comm,&eis->shell);
166: MatSetSizes(eis->shell,m,N,M,N);
167: MatSetType(eis->shell,MATSHELL);
168: MatShellSetContext(eis->shell,(void*)pc);
169: PetscLogObjectParent(pc,eis->shell);
170: MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);
171: }
172: if (!eis->usediag) return(0);
173: if (!pc->setupcalled) {
174: MatGetVecs(pc->pmat,&eis->diag,0);
175: PetscLogObjectParent(pc,eis->diag);
176: }
177: MatGetDiagonal(pc->pmat,eis->diag);
178: return(0);
179: }
181: /* --------------------------------------------------------------------*/
186: PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
187: {
188: PC_Eisenstat *eis;
191: if (omega >= 2.0 || omega <= 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
192: eis = (PC_Eisenstat*)pc->data;
193: eis->omega = omega;
194: return(0);
195: }
201: PetscErrorCode PCEisenstatNoDiagonalScaling_Eisenstat(PC pc)
202: {
203: PC_Eisenstat *eis;
206: eis = (PC_Eisenstat*)pc->data;
207: eis->usediag = PETSC_FALSE;
208: return(0);
209: }
214: /*@
215: PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
216: to use with Eisenstat's trick (where omega = 1.0 by default).
218: Collective on PC
220: Input Parameters:
221: + pc - the preconditioner context
222: - omega - relaxation coefficient (0 < omega < 2)
224: Options Database Key:
225: . -pc_eisenstat_omega <omega> - Sets omega
227: Notes:
228: The Eisenstat trick implementation of SSOR requires about 50% of the
229: usual amount of floating point operations used for SSOR + Krylov method;
230: however, the preconditioned problem must be solved with both left
231: and right preconditioning.
233: To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
234: which can be chosen with the database options
235: $ -pc_type sor -pc_sor_symmetric
237: Level: intermediate
239: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega
241: .seealso: PCSORSetOmega()
242: @*/
243: PetscErrorCode PCEisenstatSetOmega(PC pc,PetscReal omega)
244: {
245: PetscErrorCode ierr,(*f)(PC,PetscReal);
249: PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatSetOmega_C",(void (**)(void))&f);
250: if (f) {
251: (*f)(pc,omega);
252: }
253: return(0);
254: }
258: /*@
259: PCEisenstatNoDiagonalScaling - Causes the Eisenstat preconditioner
260: not to do additional diagonal preconditioning. For matrices with a constant
261: along the diagonal, this may save a small amount of work.
263: Collective on PC
265: Input Parameter:
266: . pc - the preconditioner context
268: Options Database Key:
269: . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()
271: Level: intermediate
273: Note:
274: If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
275: likley want to use this routine since it will save you some unneeded flops.
277: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR
279: .seealso: PCEisenstatSetOmega()
280: @*/
281: PetscErrorCode PCEisenstatNoDiagonalScaling(PC pc)
282: {
283: PetscErrorCode ierr,(*f)(PC);
287: PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",(void (**)(void))&f);
288: if (f) {
289: (*f)(pc);
290: }
291: return(0);
292: }
294: /* --------------------------------------------------------------------*/
296: /*MC
297: PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
298: preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.
300: Options Database Keys:
301: + -pc_eisenstat_omega <omega> - Sets omega
302: - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()
304: Level: beginner
306: Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick
308: Notes: Only implemented for the SeqAIJ matrix format.
309: Not a true parallel SOR, in parallel this implementation corresponds to block
310: Jacobi with SOR on each block.
312: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
313: PCEisenstatNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
314: M*/
319: PetscErrorCode PCCreate_Eisenstat(PC pc)
320: {
322: PC_Eisenstat *eis;
325: PetscNew(PC_Eisenstat,&eis);
326: PetscLogObjectMemory(pc,sizeof(PC_Eisenstat));
328: pc->ops->apply = PCApply_Eisenstat;
329: pc->ops->presolve = PCPreSolve_Eisenstat;
330: pc->ops->postsolve = PCPostSolve_Eisenstat;
331: pc->ops->applyrichardson = 0;
332: pc->ops->setfromoptions = PCSetFromOptions_Eisenstat;
333: pc->ops->destroy = PCDestroy_Eisenstat;
334: pc->ops->view = PCView_Eisenstat;
335: pc->ops->setup = PCSetUp_Eisenstat;
337: pc->data = (void*)eis;
338: eis->omega = 1.0;
339: eis->b = 0;
340: eis->diag = 0;
341: eis->usediag = PETSC_TRUE;
343: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatSetOmega_C","PCEisenstatSetOmega_Eisenstat",
344: PCEisenstatSetOmega_Eisenstat);
345: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",
346: "PCEisenstatNoDiagonalScaling_Eisenstat",
347: PCEisenstatNoDiagonalScaling_Eisenstat);
348: return(0);
349: }