Actual source code: snesmfjdef.c
1: /*$Id: snesmfjdef.c,v 1.31 2001/08/30 21:37:18 balay Exp $*/
2: /*
3: Implements the default PETSc approach for computing the h
4: parameter used with the finite difference based matrix-free
5: Jacobian-vector products.
7: To make your own: clone this file and modify for your needs.
9: Mandatory functions:
10: -------------------
11: MatSNESMFCompute_ - for a given point and direction computes h
13: MatSNESMFCreate_ - fills in the MatSNESMF data structure
14: for this particular implementation
16:
17: Optional functions:
18: -------------------
19: MatSNESMFView_ - prints information about the parameters being used.
20: This is called when SNESView() or -snes_view is used.
22: MatSNESMFSetFromOptions_ - checks the options database for options that
23: apply to this method.
25: MatSNESMFDestroy_ - frees any space allocated by the routines above
27: */
29: /*
30: This include file defines the data structure MatSNESMF that
31: includes information about the computation of h. It is shared by
32: all implementations that people provide
33: */
34: #include src/mat/matimpl.h
35: #include src/snes/mf/snesmfj.h
37: /*
38: The default method has one parameter that is used to
39: "cutoff" very small values. This is stored in a data structure
40: that is only visible to this file. If your method has no parameters
41: it can omit this, if it has several simply reorganize the data structure.
42: The data structure is "hung-off" the MatSNESMF data structure in
43: the void *hctx; field.
44: */
45: typedef struct {
46: PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */
47: } MatSNESMFDefault;
49: #undef __FUNCT__
51: /*
52: MatSNESMFCompute_Default - Standard PETSc code for computing the
53: differencing paramter (h) for use with matrix-free finite differences.
55: Input Parameters:
56: + ctx - the matrix free context
57: . U - the location at which you want the Jacobian
58: - a - the direction you want the derivative
60: Output Parameter:
61: . h - the scale computed
63: */
64: static int MatSNESMFCompute_Default(MatSNESMFCtx ctx,Vec U,Vec a,PetscScalar *h)
65: {
66: MatSNESMFDefault *hctx = (MatSNESMFDefault*)ctx->hctx;
67: PetscReal nrm,sum,umin = hctx->umin;
68: PetscScalar dot;
69: int ierr;
72: if (!(ctx->count % ctx->recomputeperiod)) {
73: /*
74: This algorithm requires 2 norms and 1 inner product. Rather than
75: use directly the VecNorm() and VecDot() routines (and thus have
76: three separate collective operations, we use the VecxxxBegin/End() routines
77: */
78: VecDotBegin(U,a,&dot);
79: VecNormBegin(a,NORM_1,&sum);
80: VecNormBegin(a,NORM_2,&nrm);
81: VecDotEnd(U,a,&dot);
82: VecNormEnd(a,NORM_1,&sum);
83: VecNormEnd(a,NORM_2,&nrm);
85: /*
86: Safeguard for step sizes that are "too small"
87: */
88: if (!sum) {dot = 1.0; nrm = 1.0;}
89: #if defined(PETSC_USE_COMPLEX)
90: else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
91: else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
92: #else
93: else if (dot < umin*sum && dot >= 0.0) dot = umin*sum;
94: else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum;
95: #endif
96: *h = ctx->error_rel*dot/(nrm*nrm);
97: } else {
98: *h = ctx->currenth;
99: }
100: if (*h != *h) SETERRQ3(1,"Differencing parameter is not a number sum = %g dot = %g norm = %g",sum,PetscRealPart(dot),nrm);
101: ctx->count++;
102: return(0);
103: }
105: #undef __FUNCT__
107: /*
108: MatSNESMFView_Default - Prints information about this particular
109: method for computing h. Note that this does not print the general
110: information about the matrix-free method, as such info is printed
111: by the calling routine.
113: Input Parameters:
114: + ctx - the matrix free context
115: - viewer - the PETSc viewer
116: */
117: static int MatSNESMFView_Default(MatSNESMFCtx ctx,PetscViewer viewer)
118: {
119: MatSNESMFDefault *hctx = (MatSNESMFDefault *)ctx->hctx;
120: int ierr;
121: PetscTruth isascii;
124: /*
125: Currently this only handles the ascii file viewers, others
126: could be added, but for this type of object other viewers
127: make less sense
128: */
129: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
130: if (isascii) {
131: PetscViewerASCIIPrintf(viewer," umin=%g (minimum iterate parameter)n",hctx->umin);
132: } else {
133: SETERRQ1(1,"Viewer type %s not supported for this SNES matrix free matrix",((PetscObject)viewer)->type_name);
134: }
135: return(0);
136: }
138: #undef __FUNCT__
140: /*
141: MatSNESMFSetFromOptions_Default - Looks in the options database for
142: any options appropriate for this method.
144: Input Parameter:
145: . ctx - the matrix free context
147: */
148: static int MatSNESMFSetFromOptions_Default(MatSNESMFCtx ctx)
149: {
150: int ierr;
151: MatSNESMFDefault *hctx = (MatSNESMFDefault*)ctx->hctx;
154: PetscOptionsHead("Default matrix free parameters");
155: PetscOptionsReal("-snes_mf_umin","umin","MatSNESMFDefaultSetUmin",hctx->umin,&hctx->umin,0);
156: PetscOptionsTail();
157: return(0);
158: }
160: #undef __FUNCT__
162: /*
163: MatSNESMFDestroy_Default - Frees the space allocated by
164: MatSNESMFCreate_Default().
166: Input Parameter:
167: . ctx - the matrix free context
169: Notes:
170: Does not free the ctx, that is handled by the calling routine
171: */
172: static int MatSNESMFDestroy_Default(MatSNESMFCtx ctx)
173: {
177: PetscFree(ctx->hctx);
178: return(0);
179: }
181: EXTERN_C_BEGIN
182: #undef __FUNCT__
184: /*
185: The following two routines use the PetscObjectCompose() and PetscObjectQuery()
186: mechanism to allow the user to change the Umin parameter used in this method.
187: */
188: int MatSNESMFDefaultSetUmin_Private(Mat mat,PetscReal umin)
189: {
190: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
191: MatSNESMFDefault *hctx;
194: if (!ctx) {
195: SETERRQ(1,"MatSNESMFDefaultSetUmin() attached to non-shell matrix");
196: }
197: hctx = (MatSNESMFDefault*)ctx->hctx;
198: hctx->umin = umin;
199: return(0);
200: }
201: EXTERN_C_END
203: #undef __FUNCT__
205: /*@
206: MatSNESMFDefaultSetUmin - Sets the "umin" parameter used by the default
207: PETSc routine for computing the differencing parameter, h, which is used
208: for matrix-free Jacobian-vector products.
210: Input Parameters:
211: + A - the matrix created with MatCreateSNESMF()
212: - umin - the parameter
214: Level: advanced
216: Notes:
217: See the manual page for MatCreateSNESMF() for a complete description of the
218: algorithm used to compute h.
220: .seealso: MatSNESMFSetFunctionError(), MatCreateSNESMF()
222: @*/
223: int MatSNESMFDefaultSetUmin(Mat A,PetscReal umin)
224: {
225: int ierr,(*f)(Mat,PetscReal);
229: PetscObjectQueryFunction((PetscObject)A,"MatSNESMFDefaultSetUmin_C",(void (**)(void))&f);
230: if (f) {
231: (*f)(A,umin);
232: }
233: return(0);
234: }
236: EXTERN_C_BEGIN
237: #undef __FUNCT__
239: /*
240: MatSNESMFCreate_Default - Standard PETSc code for
241: computing h with matrix-free finite differences.
243: Input Parameter:
244: . ctx - the matrix free context created by MatSNESMFCreate()
246: */
247: int MatSNESMFCreate_Default(MatSNESMFCtx ctx)
248: {
249: MatSNESMFDefault *hctx;
250: int ierr;
254: /* allocate my own private data structure */
255: ierr = PetscNew(MatSNESMFDefault,&hctx);
256: ctx->hctx = (void*)hctx;
257: /* set a default for my parameter */
258: hctx->umin = 1.e-6;
260: /* set the functions I am providing */
261: ctx->ops->compute = MatSNESMFCompute_Default;
262: ctx->ops->destroy = MatSNESMFDestroy_Default;
263: ctx->ops->view = MatSNESMFView_Default;
264: ctx->ops->setfromoptions = MatSNESMFSetFromOptions_Default;
266: PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFDefaultSetUmin_C",
267: "MatSNESMFDefaultSetUmin_Private",
268: MatSNESMFDefaultSetUmin_Private);
269: return(0);
270: }
271: EXTERN_C_END