Actual source code: snesdnest.c


  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: PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *);

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

 16:   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};

 18:   /* System generated locals */
 19:   PetscInt i__1;
 20:   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 output 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 a PetscInt 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 a PetscInt 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__ * 2);
123:   est3 = (fval[mh + 3] - fval[mh] * 2 + fval[mh - 3]) / (*h__ * 3) / (*h__ * 3);
124:   est4 = (est1 + est2 + est3) / 3;
125:   /* Computing MAX */
126:   /* Computing PETSCMAX */
127:   d__3 = PetscMax(est1, est2);
128:   /* Computing MIN */
129:   d__4 = PetscMin(est1, est2);
130:   d__1 = PetscMax(d__3, est3) - est4;
131:   d__2 = est4 - PetscMin(d__4, est3);
132:   err2 = PetscMax(d__1, d__2);
133:   /*      write (2,123) est1, est2, est3 */
134:   /* 123  format ('Second derivative estimates', 3d12.2) */
135:   if (err2 <= PetscAbsScalar(est4) * .1) *fder2 = est4;
136:   else if (err2 < PetscAbsScalar(est4)) *fder2 = est3;
137:   else *fder2 = 0.;

139:   /*     Compute the range of function values. */
140:   f_min = fval[1];
141:   f_max = fval[1];
142:   i__1  = *nf;
143:   for (i__ = 2; i__ <= i__1; ++i__) {
144:     /* Computing MIN */
145:     d__1  = f_min;
146:     d__2  = fval[i__];
147:     f_min = PetscMin(d__1, d__2);

149:     /* Computing MAX */
150:     d__1  = f_max;
151:     d__2  = fval[i__];
152:     f_max = PetscMax(d__1, d__2);
153:   }
154:   /*     Construct the difference table. */
155:   dnoise = FALSE_;
156:   for (j = 1; j <= 6; ++j) {
157:     dsgn[j - 1]   = FALSE_;
158:     cancel[j - 1] = FALSE_;
159:     scale         = 0.;
160:     i__1          = *nf - j;
161:     for (i__ = 1; i__ <= i__1; ++i__) {
162:       fval[i__] = fval[i__ + 1] - fval[i__];
163:       if (fval[i__] == 0.) cancel[j - 1] = TRUE_;

165:       /* Computing MAX */
166:       d__1  = fval[i__];
167:       d__2  = scale;
168:       d__3  = PetscAbsScalar(d__1);
169:       scale = PetscMax(d__2, d__3);
170:     }

172:     /*        Compute the estimates for the noise level. */
173:     if (scale == 0.) 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__];
190:       d__2 = fval[i__ + 1];
191:       /* Computing MAX */
192:       d__3 = fval[i__];
193:       d__4 = fval[i__ + 1];
194:       if (PetscMin(d__1, d__2) < 0. && PetscMax(d__3, d__4) > 0.) dsgn[j - 1] = TRUE_;
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) *info = 2;
202:   else /* if (complicated condition) */ {
203:     /* Computing MIN */
204:     d__1 = PetscAbsScalar(f_max);
205:     d__2 = PetscAbsScalar(f_min);
206:     if (f_max - f_min > PetscMin(d__1, d__2) * .1) *info = 3;
207:   }
208:   if (*info != 0) return 0;

210:   /*     Determine the noise level. */
211:   /* Computing MIN */
212:   d__1 = PetscMin(eps[4], eps[5]);
213:   emin = PetscMin(d__1, eps[6]);

215:   /* Computing MAX */
216:   d__1 = PetscMax(eps[4], eps[5]);
217:   emax = PetscMax(d__1, eps[6]);

219:   if (emax <= emin * 4 && dnoise) {
220:     *fnoise = (eps[4] + eps[5] + eps[6]) / 3;
221:     if (*fder2 != 0.) {
222:       *info = 1;
223:       *hopt = PetscSqrtScalar(*fnoise / PetscAbsScalar(*fder2)) * 1.68;
224:     } else {
225:       *info = 4;
226:       *hopt = *h__ * 10;
227:     }
228:     return 0;
229:   }

231:   /* Computing MIN */
232:   d__1 = PetscMin(eps[3], eps[4]);
233:   emin = PetscMin(d__1, eps[5]);

235:   /* Computing MAX */
236:   d__1 = PetscMax(eps[3], eps[4]);
237:   emax = PetscMax(d__1, eps[5]);

239:   if (emax <= emin * 4 && dnoise) {
240:     *fnoise = (eps[3] + eps[4] + eps[5]) / 3;
241:     if (*fder2 != 0.) {
242:       *info = 1;
243:       *hopt = PetscSqrtScalar(*fnoise / PetscAbsScalar(*fder2)) * 1.68;
244:     } else {
245:       *info = 4;
246:       *hopt = *h__ * 10;
247:     }
248:     return 0;
249:   }
250:   /*     Noise not detected; decide if h is too small or too large. */
251:   if (!cancel[3]) {
252:     if (dsgn[3]) *info = 2;
253:     else *info = 3;
254:     return 0;
255:   }
256:   if (!cancel[2]) {
257:     if (dsgn[2]) *info = 2;
258:     else *info = 3;
259:     return 0;
260:   }
261:   /*     If there is cancelllation on the third and fourth column */
262:   /*     then h is too small */
263:   *info = 2;
264:   return 0;
265:   /*      if (cancel .or. dsgn(3)) then */
266:   /*         info = 2 */
267:   /*      else */
268:   /*         info = 3 */
269:   /*      end if */
270: } /* dnest_ */