1 /* fnoise/snesdnest.F -- translated by f2c (version 20020314).
2 */
3 #include <petscsys.h>
4 #define FALSE_ 0
5 #define TRUE_ 1
6
7 /* Noise estimation routine, written by Jorge More'. Details are below. */
8
9 PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *);
10
SNESNoise_dnest_(PetscInt * nf,double * fval,double * h__,double * fnoise,double * fder2,double * hopt,PetscInt * info,double * eps)11 PetscErrorCode SNESNoise_dnest_(PetscInt *nf, double *fval, double *h__, double *fnoise, double *fder2, double *hopt, PetscInt *info, double *eps)
12 {
13 /* Initialized data */
14
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};
16
17 /* System generated locals */
18 PetscInt i__1;
19 double d__1, d__2, d__3, d__4;
20
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;
30
31 /* ********** */
32
33 /* Subroutine dnest */
34
35 /* This subroutine estimates the noise in a function */
36 /* and provides estimates of the optimal difference parameter */
37 /* for a forward-difference approximation. */
38
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 */
42
43 /* f(x-2*h), f(x-h), f(x), f(x+h), f(x+2*h), */
44
45 /* in the array fval. The use of nf = 7 function evaluations is */
46 /* recommended. */
47
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. */
52
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. */
58
59 /* The subroutine statement is */
60
61 /* subroutine dnest(nf,fval,h,hopt,fnoise,info,eps) */
62
63 /* where */
64
65 /* nf is a PetscInt variable. */
66 /* On entry nf is the number of function values. */
67 /* On exit nf is unchanged. */
68
69 /* f is a double precision array of dimension nf. */
70 /* On entry f contains the function values. */
71 /* On exit f is overwritten. */
72
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. */
76
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. */
81
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. */
86
87 /* info is a PetscInt variable. */
88 /* On entry info need not be specified. */
89 /* On exit info is set as follows: */
90
91 /* info = 1 Noise has been detected. */
92
93 /* info = 2 Noise has not been detected; h is too small. */
94 /* Try 100*h for the next value of h. */
95
96 /* info = 3 Noise has not been detected; h is too large. */
97 /* Try h/100 for the next value of h. */
98
99 /* info = 4 Noise has been detected but the estimate of hopt */
100 /* is not reliable; h is too small. */
101
102 /* eps is a double precision work array of dimension nf. */
103
104 /* MINPACK-2 Project. April 1997. */
105 /* Argonne National Laboratory. */
106 /* Jorge J. More'. */
107
108 /* ********** */
109 /* Parameter adjustments */
110 --eps;
111 --fval;
112
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.;
137
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);
147
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_;
163
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 }
170
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);
208
209 /* Determine the noise level. */
210 /* Computing MIN */
211 d__1 = PetscMin(eps[4], eps[5]);
212 emin = PetscMin(d__1, eps[6]);
213
214 /* Computing MAX */
215 d__1 = PetscMax(eps[4], eps[5]);
216 emax = PetscMax(d__1, eps[6]);
217
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 }
229
230 /* Computing MIN */
231 d__1 = PetscMin(eps[3], eps[4]);
232 emin = PetscMin(d__1, eps[5]);
233
234 /* Computing MAX */
235 d__1 = PetscMax(eps[3], eps[4]);
236 emax = PetscMax(d__1, eps[5]);
237
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_ */
270