Actual source code: snesdnest.c
1: /* fnoise/snesdnest.F -- translated by f2c (version 20020314).
2: */
3: #include <petscsys.h>
4: #define FALSE_ 0
5: #define TRUE_ 1
7: /* Noise estimation routine, written by Jorge More'. Details are below. */
9: PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *);
11: PetscErrorCode SNESNoise_dnest_(PetscInt *nf, double *fval, double *h__, double *fnoise, double *fder2, double *hopt, PetscInt *info, double *eps)
12: {
13: /* Initialized data */
15: static double const__[15] = {.71, .41, .23, .12, .063, .033, .018, .0089, .0046, .0024, .0012, 6.1e-4, 3.1e-4, 1.6e-4, 8e-5};
17: /* System generated locals */
18: PetscInt i__1;
19: double d__1, d__2, d__3, d__4;
21: /* Local variables */
22: static double emin, emax;
23: static PetscInt dsgn[6];
24: static double f_max, f_min, stdv;
25: static PetscInt i__, j;
26: static double scale;
27: static PetscInt mh;
28: static PetscInt cancel[6], dnoise;
29: static double err2, est1, est2, est3, est4;
31: /* ********** */
33: /* Subroutine dnest */
35: /* This subroutine estimates the noise in a function */
36: /* and provides estimates of the optimal difference parameter */
37: /* for a forward-difference approximation. */
39: /* The user must provide a difference parameter h, and the */
40: /* function value at nf points centered around the current point. */
41: /* For example, if nf = 7, the user must provide */
43: /* f(x-2*h), f(x-h), f(x), f(x+h), f(x+2*h), */
45: /* in the array fval. The use of nf = 7 function evaluations is */
46: /* recommended. */
48: /* The noise in the function is roughly defined as the variance in */
49: /* the computed value of the function. The noise in the function */
50: /* provides valuable information. For example, function values */
51: /* smaller than the noise should be considered to be zero. */
53: /* This subroutine requires an initial estimate for h. Under estimates */
54: /* are usually preferred. If noise is not detected, the user should */
55: /* increase or decrease h according to the output value of info. */
56: /* In most cases, the subroutine detects noise with the initial */
57: /* value of h. */
59: /* The subroutine statement is */
61: /* subroutine dnest(nf,fval,h,hopt,fnoise,info,eps) */
63: /* where */
65: /* nf is a PetscInt variable. */
66: /* On entry nf is the number of function values. */
67: /* On exit nf is unchanged. */
69: /* f is a double precision array of dimension nf. */
70: /* On entry f contains the function values. */
71: /* On exit f is overwritten. */
73: /* h is a double precision variable. */
74: /* On entry h is an estimate of the optimal difference parameter. */
75: /* On exit h is unchanged. */
77: /* fnoise is a double precision variable. */
78: /* On entry fnoise need not be specified. */
79: /* On exit fnoise is set to an estimate of the function noise */
80: /* if noise is detected; otherwise fnoise is set to zero. */
82: /* hopt is a double precision variable. */
83: /* On entry hopt need not be specified. */
84: /* On exit hopt is set to an estimate of the optimal difference */
85: /* parameter if noise is detected; otherwise hopt is set to zero. */
87: /* info is a PetscInt variable. */
88: /* On entry info need not be specified. */
89: /* On exit info is set as follows: */
91: /* info = 1 Noise has been detected. */
93: /* info = 2 Noise has not been detected; h is too small. */
94: /* Try 100*h for the next value of h. */
96: /* info = 3 Noise has not been detected; h is too large. */
97: /* Try h/100 for the next value of h. */
99: /* info = 4 Noise has been detected but the estimate of hopt */
100: /* is not reliable; h is too small. */
102: /* eps is a double precision work array of dimension nf. */
104: /* MINPACK-2 Project. April 1997. */
105: /* Argonne National Laboratory. */
106: /* Jorge J. More'. */
108: /* ********** */
109: /* Parameter adjustments */
110: --eps;
111: --fval;
113: /* Function Body */
114: *fnoise = 0.;
115: *fder2 = 0.;
116: *hopt = 0.;
117: /* Compute an estimate of the second derivative and */
118: /* determine a bound on the error. */
119: mh = (*nf + 1) / 2;
120: est1 = (fval[mh + 1] - fval[mh] * 2 + fval[mh - 1]) / *h__ / *h__;
121: est2 = (fval[mh + 2] - fval[mh] * 2 + fval[mh - 2]) / (*h__ * 2) / (*h__ * 2);
122: est3 = (fval[mh + 3] - fval[mh] * 2 + fval[mh - 3]) / (*h__ * 3) / (*h__ * 3);
123: est4 = (est1 + est2 + est3) / 3;
124: /* Computing MAX */
125: /* Computing PETSCMAX */
126: d__3 = PetscMax(est1, est2);
127: /* Computing MIN */
128: d__4 = PetscMin(est1, est2);
129: d__1 = PetscMax(d__3, est3) - est4;
130: d__2 = est4 - PetscMin(d__4, est3);
131: err2 = PetscMax(d__1, d__2);
132: /* write (2,123) est1, est2, est3 */
133: /* 123 format ('Second derivative estimates', 3d12.2) */
134: if (err2 <= PetscAbsScalar(est4) * .1) *fder2 = est4;
135: else if (err2 < PetscAbsScalar(est4)) *fder2 = est3;
136: else *fder2 = 0.;
138: /* Compute the range of function values. */
139: f_min = fval[1];
140: f_max = fval[1];
141: i__1 = *nf;
142: for (i__ = 2; i__ <= i__1; ++i__) {
143: /* Computing MIN */
144: d__1 = f_min;
145: d__2 = fval[i__];
146: f_min = PetscMin(d__1, d__2);
148: /* Computing MAX */
149: d__1 = f_max;
150: d__2 = fval[i__];
151: f_max = PetscMax(d__1, d__2);
152: }
153: /* Construct the difference table. */
154: dnoise = FALSE_;
155: for (j = 1; j <= 6; ++j) {
156: dsgn[j - 1] = FALSE_;
157: cancel[j - 1] = FALSE_;
158: scale = 0.;
159: i__1 = *nf - j;
160: for (i__ = 1; i__ <= i__1; ++i__) {
161: fval[i__] = fval[i__ + 1] - fval[i__];
162: if (fval[i__] == 0.) cancel[j - 1] = TRUE_;
164: /* Computing MAX */
165: d__1 = fval[i__];
166: d__2 = scale;
167: d__3 = PetscAbsScalar(d__1);
168: scale = PetscMax(d__2, d__3);
169: }
171: /* Compute the estimates for the noise level. */
172: if (scale == 0.) stdv = 0.;
173: else {
174: stdv = 0.;
175: i__1 = *nf - j;
176: for (i__ = 1; i__ <= i__1; ++i__) {
177: /* Computing 2nd power */
178: d__1 = fval[i__] / scale;
179: stdv += d__1 * d__1;
180: }
181: stdv = scale * PetscSqrtScalar(stdv / (*nf - j));
182: }
183: eps[j] = const__[j - 1] * stdv;
184: /* Determine differences in sign. */
185: i__1 = *nf - j - 1;
186: for (i__ = 1; i__ <= i__1; ++i__) {
187: /* Computing MIN */
188: d__1 = fval[i__];
189: d__2 = fval[i__ + 1];
190: /* Computing MAX */
191: d__3 = fval[i__];
192: d__4 = fval[i__ + 1];
193: if (PetscMin(d__1, d__2) < 0. && PetscMax(d__3, d__4) > 0.) dsgn[j - 1] = TRUE_;
194: }
195: }
196: /* First requirement for detection of noise. */
197: dnoise = dsgn[3];
198: /* Check for h too small or too large. */
199: *info = 0;
200: if (f_max == f_min) *info = 2;
201: else /* if (complicated condition) */ {
202: /* Computing MIN */
203: d__1 = PetscAbsScalar(f_max);
204: d__2 = PetscAbsScalar(f_min);
205: if (f_max - f_min > PetscMin(d__1, d__2) * .1) *info = 3;
206: }
207: if (*info != 0) PetscFunctionReturn(PETSC_SUCCESS);
209: /* Determine the noise level. */
210: /* Computing MIN */
211: d__1 = PetscMin(eps[4], eps[5]);
212: emin = PetscMin(d__1, eps[6]);
214: /* Computing MAX */
215: d__1 = PetscMax(eps[4], eps[5]);
216: emax = PetscMax(d__1, eps[6]);
218: if (emax <= emin * 4 && dnoise) {
219: *fnoise = (eps[4] + eps[5] + eps[6]) / 3;
220: if (*fder2 != 0.) {
221: *info = 1;
222: *hopt = PetscSqrtScalar(*fnoise / PetscAbsScalar(*fder2)) * 1.68;
223: } else {
224: *info = 4;
225: *hopt = *h__ * 10;
226: }
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /* Computing MIN */
231: d__1 = PetscMin(eps[3], eps[4]);
232: emin = PetscMin(d__1, eps[5]);
234: /* Computing MAX */
235: d__1 = PetscMax(eps[3], eps[4]);
236: emax = PetscMax(d__1, eps[5]);
238: if (emax <= emin * 4 && dnoise) {
239: *fnoise = (eps[3] + eps[4] + eps[5]) / 3;
240: if (*fder2 != 0.) {
241: *info = 1;
242: *hopt = PetscSqrtScalar(*fnoise / PetscAbsScalar(*fder2)) * 1.68;
243: } else {
244: *info = 4;
245: *hopt = *h__ * 10;
246: }
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
249: /* Noise not detected; decide if h is too small or too large. */
250: if (!cancel[3]) {
251: if (dsgn[3]) *info = 2;
252: else *info = 3;
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
255: if (!cancel[2]) {
256: if (dsgn[2]) *info = 2;
257: else *info = 3;
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
260: /* If there is cancellation on the third and fourth column */
261: /* then h is too small */
262: *info = 2;
263: PetscFunctionReturn(PETSC_SUCCESS);
264: /* if (cancel .or. dsgn(3)) then */
265: /* info = 2 */
266: /* else */
267: /* info = 3 */
268: /* end if */
269: } /* dnest_ */