2: /* fnoise/snesdnest.F -- translated by f2c (version 20020314).
  3: */
  4: #include <petscsys.h>
  5: #define FALSE_ 0
  6: #define TRUE_ 1

  8: /*  Noise estimation routine, written by Jorge More'.  Details are below. */

 10: /* Subroutine */ PetscErrorCode SNESNoise_dnest_(PetscInt *nf, double *fval,double *h__,double *fnoise, double *fder2, double *hopt, PetscInt *info, double *eps)
 11: {
 12:     /* Initialized data */

 14:     static double const__[15] = { .71,.41,.23,.12,.063,.033,.018,.0089,
 15:             .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;


 22:     /* Local variables */
 23:     static double emin, emax;
 24:     static PetscInt dsgn[6];
 25:     static double f_max, f_min, stdv;
 26:     static PetscInt i__, j;
 27:     static double scale;
 28:     static PetscInt mh;
 29:     static PetscInt cancel[6], dnoise;
 30:     static double err2, est1, est2, est3, est4;

 32: /*     ********** */

 34: /*     Subroutine dnest */

 36: /*     This subroutine estimates the noise in a function */
 37: /*     and provides estimates of the optimal difference parameter */
 38: /*     for a forward-difference approximation. */

 40: /*     The user must provide a difference parameter h, and the */
 41: /*     function value at nf points centered around the current point. */
 42: /*     For example, if nf = 7, the user must provide */

 44: /*        f(x-2*h), f(x-h), f(x), f(x+h),  f(x+2*h), */

 46: /*     in the array fval. The use of nf = 7 function evaluations is */
 47: /*     recommended. */

 49: /*     The noise in the function is roughly defined as the variance in */
 50: /*     the computed value of the function. The noise in the function */
 51: /*     provides valuable information. For example, function values */
 52: /*     smaller than the noise should be considered to be zero. */

 54: /*     This subroutine requires an initial estimate for h. Under estimates */
 55: /*     are usually preferred. If noise is not detected, the user should */
 56: /*     increase or decrease h according to the ouput value of info. */
 57: /*     In most cases, the subroutine detects noise with the initial */
 58: /*     value of h. */

 60: /*     The subroutine statement is */

 62: /*       subroutine dnest(nf,fval,h,hopt,fnoise,info,eps) */

 64: /*     where */

 66: /*       nf is an int variable. */
 67: /*         On entry nf is the number of function values. */
 68: /*         On exit nf is unchanged. */

 70: /*       f is a double precision array of dimension nf. */
 71: /*         On entry f contains the function values. */
 72: /*         On exit f is overwritten. */

 74: /*       h is a double precision variable. */
 75: /*         On entry h is an estimate of the optimal difference parameter. */
 76: /*         On exit h is unchanged. */

 78: /*       fnoise is a double precision variable. */
 79: /*         On entry fnoise need not be specified. */
 80: /*         On exit fnoise is set to an estimate of the function noise */
 81: /*            if noise is detected; otherwise fnoise is set to zero. */

 83: /*       hopt is a double precision variable. */
 84: /*         On entry hopt need not be specified. */
 85: /*         On exit hopt is set to an estimate of the optimal difference */
 86: /*            parameter if noise is detected; otherwise hopt is set to zero. */

 88: /*       info is an int variable. */
 89: /*         On entry info need not be specified. */
 90: /*         On exit info is set as follows: */

 92: /*            info = 1  Noise has been detected. */

 94: /*            info = 2  Noise has not been detected; h is too small. */
 95: /*                      Try 100*h for the next value of h. */

 97: /*            info = 3  Noise has not been detected; h is too large. */
 98: /*                      Try h/100 for the next value of h. */

100: /*            info = 4  Noise has been detected but the estimate of hopt */
101: /*                      is not reliable; h is too small. */

103: /*       eps is a double precision work array of dimension nf. */

105: /*     MINPACK-2 Project. April 1997. */
106: /*     Argonne National Laboratory. */
107: /*     Jorge J. More'. */

109: /*     ********** */
110:     /* Parameter adjustments */
111:     --eps;
112:     --fval;

114:     /* Function Body */
115:     *fnoise = 0.;
116:     *fder2 = 0.;
117:     *hopt = 0.;
118: /*     Compute an estimate of the second derivative and */
119: /*     determine a bound on the error. */
120:     mh = (*nf + 1) / 2;
121:     est1 = (fval[mh + 1] - fval[mh] * 2 + fval[mh - 1]) / *h__ / *h__;
122:     est2 = (fval[mh + 2] - fval[mh] * 2 + fval[mh - 2]) / (*h__ * 2) / (*h__ *
123:              2);
124:     est3 = (fval[mh + 3] - fval[mh] * 2 + fval[mh - 3]) / (*h__ * 3) / (*h__ *
125:              3);
126:     est4 = (est1 + est2 + est3) / 3;
127: /* Computing MAX */
128: /* Computing PETSCMAX */
129:     d__3 = PetscMax(est1,est2);
130: /* Computing MIN */
131:     d__4 = PetscMin(est1,est2);
132:     d__1 = PetscMax(d__3,est3) - est4, d__2 = est4 - PetscMin(d__4,est3);
133:     err2 = PetscMax(d__1,d__2);
134: /*      write (2,123) est1, est2, est3 */
135: /* 123  format ('Second derivative estimates', 3d12.2) */
136:     if (err2 <= PetscAbsScalar(est4) * .1) {
137:         *fder2 = est4;
138:     } else if (err2 < PetscAbsScalar(est4)) {
139:         *fder2 = est3;
140:     } else {
141:         *fder2 = 0.;
142:     }
143: /*     Compute the range of function values. */
144:     f_min = fval[1];
145:     f_max = fval[1];
146:     i__1 = *nf;
147:     for (i__ = 2; i__ <= i__1; ++i__) {
148: /* Computing MIN */
149:         d__1 = f_min, d__2 = fval[i__];
150:         f_min = PetscMin(d__1,d__2);
151: /* Computing MAX */
152:         d__1 = f_max, d__2 = fval[i__];
153:         f_max = PetscMax(d__1,d__2);
154:     }
155: /*     Construct the difference table. */
156:     dnoise = FALSE_;
157:     for (j = 1; j <= 6; ++j) {
158:         dsgn[j - 1] = FALSE_;
159:         cancel[j - 1] = FALSE_;
160:         scale = 0.;
161:         i__1 = *nf - j;
162:         for (i__ = 1; i__ <= i__1; ++i__) {
163:             fval[i__] = fval[i__ + 1] - fval[i__];
164:             if (fval[i__] == 0.) {
165:                 cancel[j - 1] = TRUE_;
166:             }
167: /* Computing MAX */
168:             d__2 = scale, d__3 = (d__1 = fval[i__], PetscAbsScalar(d__1));
169:             scale = PetscMax(d__2,d__3);
170:         }
171: /*        Compute the estimates for the noise level. */
172:         if (scale == 0.) {
173:             stdv = 0.;
174:         } else {
175:             stdv = 0.;
176:             i__1 = *nf - j;
177:             for (i__ = 1; i__ <= i__1; ++i__) {
178: /* Computing 2nd power */
179:                 d__1 = fval[i__] / scale;
180:                 stdv += d__1 * d__1;
181:             }
182:             stdv = scale * PetscSqrtScalar(stdv / (*nf - j));
183:         }
184:         eps[j] = const__[j - 1] * stdv;
185: /*        Determine differences in sign. */
186:         i__1 = *nf - j - 1;
187:         for (i__ = 1; i__ <= i__1; ++i__) {
188: /* Computing MIN */
189:             d__1 = fval[i__], d__2 = fval[i__ + 1];
190: /* Computing MAX */
191:             d__3 = fval[i__], d__4 = fval[i__ + 1];
192:             if (PetscMin(d__1,d__2) < 0. && PetscMax(d__3,d__4) > 0.) {
193:                 dsgn[j - 1] = TRUE_;
194:             }
195:         }
196:     }
197: /*     First requirement for detection of noise. */
198:     dnoise = dsgn[3];
199: /*     Check for h too small or too large. */
200:     *info = 0;
201:     if (f_max == f_min) {
202:         *info = 2;
203:     } else /* if(complicated condition) */ {
204: /* Computing MIN */
205:         d__1 = PetscAbsScalar(f_max), d__2 = PetscAbsScalar(f_min);
206:         if (f_max - f_min > PetscMin(d__1,d__2) * .1) {
207:             *info = 3;
208:         }
209:     }
210:     if (*info != 0) {
211:         return(0);
212:     }
213: /*     Determine the noise level. */
214: /* Computing MIN */
215:     d__1 = PetscMin(eps[4],eps[5]);
216:     emin = PetscMin(d__1,eps[6]);
217: /* Computing MAX */
218:     d__1 = PetscMax(eps[4],eps[5]);
219:     emax = PetscMax(d__1,eps[6]);
220:     if (emax <= emin * 4 && dnoise) {
221:         *fnoise = (eps[4] + eps[5] + eps[6]) / 3;
222:         if (*fder2 != 0.) {
223:             *info = 1;
224:             *hopt = PetscSqrtScalar(*fnoise / PetscAbsScalar(*fder2)) * 1.68;
225:         } else {
226:             *info = 4;
227:             *hopt = *h__ * 10;
228:         }
229:         return(0);
230:     }
231: /* Computing MIN */
232:     d__1 = PetscMin(eps[3],eps[4]);
233:     emin = PetscMin(d__1,eps[5]);
234: /* Computing MAX */
235:     d__1 = PetscMax(eps[3],eps[4]);
236:     emax = PetscMax(d__1,eps[5]);
237:     if (emax <= emin * 4 && dnoise) {
238:         *fnoise = (eps[3] + eps[4] + eps[5]) / 3;
239:         if (*fder2 != 0.) {
240:             *info = 1;
241:             *hopt = PetscSqrtScalar(*fnoise / PetscAbsScalar(*fder2)) * 1.68;
242:         } else {
243:             *info = 4;
244:             *hopt = *h__ * 10;
245:         }
246:         return(0);
247:     }
248: /*     Noise not detected; decide if h is too small or too large. */
249:     if (! cancel[3]) {
250:         if (dsgn[3]) {
251:             *info = 2;
252:         } else {
253:             *info = 3;
254:         }
255:         return(0);
256:     }
257:     if (! cancel[2]) {
258:         if (dsgn[2]) {
259:             *info = 2;
260:         } else {
261:             *info = 3;
262:         }
263:         return(0);
264:     }
265: /*     If there is cancelllation on the third and fourth column */
266: /*     then h is too small */
267:     *info = 2;
268:     return(0);
269: /*      if (cancel .or. dsgn(3)) then */
270: /*         info = 2 */
271: /*      else */
272: /*         info = 3 */
273: /*      end if */
274: } /* dnest_ */