Actual source code: wp.c
1: /*$Id: wp.c,v 1.36 2001/08/06 21:17:13 bsmith Exp $*/
2: /*
3: Implements an alternative approach for computing the differencing parameter
4: h used with the finite difference based matrix-free Jacobian. This code
5: implements the strategy of M. Pernice and H. Walker:
7: h = error_rel * sqrt(1 + ||U||) / ||a||
9: Notes:
10: 1) || U || does not change between linear iterations so can be reused
11: 2) In GMRES || a || == 1 and so does not need to ever be computed if you never
12: have a restart. Unfortunately a RESTART computes a matrix vector product
13: with ||a|| != 0 which breaks this
15: Reference: M. Pernice and H. F. Walker, "NITSOL: A Newton Iterative
16: Solver for Nonlinear Systems", SIAM J. Sci. Stat. Comput.", 1998,
17: vol 19, pp. 302--318.
19: See snesmfjdef.c for a full set of comments on the routines below.
20: */
22: /*
23: This include file defines the data structure MatSNESMF that
24: includes information about the computation of h. It is shared by
25: all implementations that people provide.
26: */
27: #include src/mat/matimpl.h
28: #include src/snes/mf/snesmfj.h
30: typedef struct {
31: PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */
32: PetscTruth computenorma,computenormU;
33: } MatSNESMFWP;
35: #undef __FUNCT__
37: /*
38: MatSNESMFCompute_WP - Standard PETSc code for
39: computing h with matrix-free finite differences.
41: Input Parameters:
42: + ctx - the matrix free context
43: . U - the location at which you want the Jacobian
44: - a - the direction you want the derivative
46: Output Parameter:
47: . h - the scale computed
49: */
50: static int MatSNESMFCompute_WP(MatSNESMFCtx ctx,Vec U,Vec a,PetscScalar *h)
51: {
52: MatSNESMFWP *hctx = (MatSNESMFWP*)ctx->hctx;
53: PetscReal normU,norma = 1.0;
54: int ierr;
58: if (!(ctx->count % ctx->recomputeperiod)) {
59: if (hctx->computenorma && (hctx->computenormU || !ctx->ncurrenth)) {
60: VecNormBegin(U,NORM_2,&normU);
61: VecNormBegin(a,NORM_2,&norma);
62: VecNormEnd(U,NORM_2,&normU);
63: VecNormEnd(a,NORM_2,&norma);
64: hctx->normUfact = sqrt(1.0+normU);
65: } else if (hctx->computenormU || !ctx->ncurrenth) {
66: VecNorm(U,NORM_2,&normU);
67: hctx->normUfact = sqrt(1.0+normU);
68: } else if (hctx->computenorma) {
69: VecNorm(a,NORM_2,&norma);
70: }
71: *h = ctx->error_rel*hctx->normUfact/norma;
72: } else {
73: *h = ctx->currenth;
74: }
75: return(0);
76: }
78: #undef __FUNCT__
80: /*
81: MatSNESMFView_WP - Prints information about this particular
82: method for computing h. Note that this does not print the general
83: information about the matrix free, that is printed by the calling
84: routine.
86: Input Parameters:
87: + ctx - the matrix free context
88: - viewer - the PETSc viewer
90: */
91: static int MatSNESMFView_WP(MatSNESMFCtx ctx,PetscViewer viewer)
92: {
93: MatSNESMFWP *hctx = (MatSNESMFWP *)ctx->hctx;
94: int ierr;
95: PetscTruth isascii;
98: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
99: if (isascii) {
100: if (hctx->computenorma){PetscViewerASCIIPrintf(viewer," Computes normAn");}
101: else { PetscViewerASCIIPrintf(viewer," Does not compute normAn");}
102: if (hctx->computenormU){ PetscViewerASCIIPrintf(viewer," Computes normUn");}
103: else { PetscViewerASCIIPrintf(viewer," Does not compute normUn");}
104: } else {
105: SETERRQ1(1,"Viewer type %s not supported for SNES matrix-free WP",((PetscObject)viewer)->type_name);
106: }
107: return(0);
108: }
110: #undef __FUNCT__
112: /*
113: MatSNESMFSetFromOptions_WP - Looks in the options database for
114: any options appropriate for this method
116: Input Parameter:
117: . ctx - the matrix free context
119: */
120: static int MatSNESMFSetFromOptions_WP(MatSNESMFCtx ctx)
121: {
122: int ierr;
123: MatSNESMFWP *hctx = (MatSNESMFWP*)ctx->hctx;
126: PetscOptionsHead("Walker-Pernice options");
127: PetscOptionsLogical("-snes_mf_compute_norma","Compute the norm of a","MatSNESMFWPSetComputeNormA",
128: hctx->computenorma,&hctx->computenorma,0);
129: PetscOptionsLogical("-snes_mf_compute_normu","Compute the norm of u","MatSNESMFWPSetComputeNormU",
130: hctx->computenorma,&hctx->computenorma,0);
131: PetscOptionsTail();
132: return(0);
133: }
135: #undef __FUNCT__
137: /*
138: MatSNESMFDestroy_WP - Frees the space allocated by
139: MatSNESMFCreate_WP().
141: Input Parameter:
142: . ctx - the matrix free context
144: Notes: does not free the ctx, that is handled by the calling routine
146: */
147: static int MatSNESMFDestroy_WP(MatSNESMFCtx ctx)
148: {
151: PetscFree(ctx->hctx);
152: return(0);
153: }
155: EXTERN_C_BEGIN
156: #undef __FUNCT__
158: int MatSNESMFWPSetComputeNormA_P(Mat mat,PetscTruth flag)
159: {
160: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
161: MatSNESMFWP *hctx;
164: hctx = (MatSNESMFWP*)ctx->hctx;
165: hctx->computenorma = flag;
166: return(0);
167: }
168: EXTERN_C_END
170: #undef __FUNCT__
172: /*@
173: MatSNESMFWPSetComputeNormA - Sets whether it computes the ||a|| used by the WP
174: PETSc routine for computing h. With GMRES since the ||a|| is always
175: one, you can save communication by setting this to false.
177: Input Parameters:
178: + A - the matrix created with MatCreateSNESMF()
179: - flag - PETSC_TRUE causes it to compute ||a||, PETSC_FALSE assumes it is 1.
181: Level: advanced
183: Notes:
184: See the manual page for MatCreateSNESMF() for a complete description of the
185: algorithm used to compute h.
187: .seealso: MatSNESMFSetFunctionError(), MatCreateSNESMF()
189: @*/
190: int MatSNESMFWPSetComputeNormA(Mat A,PetscTruth flag)
191: {
192: int ierr,(*f)(Mat,PetscTruth);
196: PetscObjectQueryFunction((PetscObject)A,"MatSNESMFWPSetComputeNormA_C",(void (**)(void))&f);
197: if (f) {
198: (*f)(A,flag);
199: }
200: return(0);
201: }
203: EXTERN_C_BEGIN
204: #undef __FUNCT__
206: int MatSNESMFWPSetComputeNormU_P(Mat mat,PetscTruth flag)
207: {
208: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
209: MatSNESMFWP *hctx;
212: hctx = (MatSNESMFWP*)ctx->hctx;
213: hctx->computenormU = flag;
214: return(0);
215: }
216: EXTERN_C_END
218: #undef __FUNCT__
220: /*@
221: MatSNESMFWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
222: PETSc routine for computing h. With any Krylov solver this need only
223: be computed during the first iteration and kept for later.
225: Input Parameters:
226: + A - the matrix created with MatCreateSNESMF()
227: - flag - PETSC_TRUE causes it to compute ||U||, PETSC_FALSE uses the previous value
229: Level: advanced
231: Notes:
232: See the manual page for MatCreateSNESMF() for a complete description of the
233: algorithm used to compute h.
235: .seealso: MatSNESMFSetFunctionError(), MatCreateSNESMF()
237: @*/
238: int MatSNESMFWPSetComputeNormU(Mat A,PetscTruth flag)
239: {
240: int ierr,(*f)(Mat,PetscTruth);
244: PetscObjectQueryFunction((PetscObject)A,"MatSNESMFWPSetComputeNormU_C",(void (**)(void))&f);
245: if (f) {
246: (*f)(A,flag);
247: }
248: return(0);
249: }
251: EXTERN_C_BEGIN
252: #undef __FUNCT__
254: /*
255: MatSNESMFCreate_WP - Standard PETSc code for
256: computing h with matrix-free finite differences.
258: Input Parameter:
259: . ctx - the matrix free context created by MatSNESMFCreate()
261: */
262: int MatSNESMFCreate_WP(MatSNESMFCtx ctx)
263: {
264: int ierr;
265: MatSNESMFWP *hctx;
269: /* allocate my own private data structure */
270: ierr = PetscNew(MatSNESMFWP,&hctx);
271: ctx->hctx = (void*)hctx;
272: hctx->computenormU = PETSC_FALSE;
273: hctx->computenorma = PETSC_TRUE;
275: /* set the functions I am providing */
276: ctx->ops->compute = MatSNESMFCompute_WP;
277: ctx->ops->destroy = MatSNESMFDestroy_WP;
278: ctx->ops->view = MatSNESMFView_WP;
279: ctx->ops->setfromoptions = MatSNESMFSetFromOptions_WP;
281: PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFWPSetComputeNormA_C",
282: "MatSNESMFWPSetComputeNormA_P",
283: MatSNESMFWPSetComputeNormA_P);
284: PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFWPSetComputeNormU_C",
285: "MatSNESMFWPSetComputeNormU_P",
286: MatSNESMFWPSetComputeNormU_P);
288: return(0);
289: }
290: EXTERN_C_END