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,
 17:                                 .0046,.0024,.0012,6.1e-4,3.1e-4,1.6e-4,8e-5 };

 19:   /* System generated locals */
 20:   PetscInt i__1;
 21:   double   d__1, d__2, d__3, d__4;

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

 33: /*     ********** */

 35: /*     Subroutine dnest */

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

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

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

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

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

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

 61: /*     The subroutine statement is */

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

 65: /*     where */

 67: /*       nf is a PetscInt variable. */
 68: /*         On entry nf is the number of function values. */
 69: /*         On exit nf is unchanged. */

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

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

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

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

 89: /*       info is a PetscInt variable. */
 90: /*         On entry info need not be specified. */
 91: /*         On exit info is set as follows: */

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

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

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

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

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

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

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

115:   /* Function Body */
116:   *fnoise = 0.;
117:   *fder2  = 0.;
118:   *hopt   = 0.;
119: /*     Compute an estimate of the second derivative and */
120: /*     determine a bound on the error. */
121:   mh   = (*nf + 1) / 2;
122:   est1 = (fval[mh + 1] - fval[mh] * 2 + fval[mh - 1]) / *h__ / *h__;
123:   est2 = (fval[mh + 2] - fval[mh] * 2 + fval[mh - 2]) / (*h__ * 2) / (*h__ * 2);
124:   est3 = (fval[mh + 3] - fval[mh] * 2 + fval[mh - 3]) / (*h__ * 3) / (*h__ * 3);
125:   est4 = (est1 + est2 + est3) / 3;
126: /* Computing MAX */
127: /* Computing PETSCMAX */
128:   d__3 = PetscMax(est1,est2);
129: /* Computing MIN */
130:   d__4 = PetscMin(est1,est2);
131:   d__1 = PetscMax(d__3,est3) - est4;
132:   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) *fder2 = est4;
137:   else if (err2 < PetscAbsScalar(est4))  *fder2 = est3;
138:   else *fder2 = 0.;

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

150:     /* Computing MAX */
151:     d__1 = f_max;
152:     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.) cancel[j - 1] = TRUE_;

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

173:     /*        Compute the estimates for the noise level. */
174:     if (scale == 0.) stdv = 0.;
175:     else {
176:       stdv = 0.;
177:       i__1 = *nf - j;
178:       for (i__ = 1; i__ <= i__1; ++i__) {
179:         /* Computing 2nd power */
180:         d__1 = fval[i__] / scale;
181:         stdv += d__1 * d__1;
182:       }
183:       stdv = scale * PetscSqrtScalar(stdv / (*nf - j));
184:     }
185:     eps[j] = const__[j - 1] * stdv;
186: /*        Determine differences in sign. */
187:     i__1 = *nf - j - 1;
188:     for (i__ = 1; i__ <= i__1; ++i__) {
189:       /* Computing MIN */
190:       d__1 = fval[i__];
191:       d__2 = fval[i__ + 1];
192:       /* Computing MAX */
193:       d__3 = fval[i__];
194:       d__4 = fval[i__ + 1];
195:       if (PetscMin(d__1,d__2) < 0. && PetscMax(d__3,d__4) > 0.) dsgn[j - 1] = TRUE_;
196:     }
197:   }
198:   /*     First requirement for detection of noise. */
199:   dnoise = dsgn[3];
200:   /*     Check for h too small or too large. */
201:   *info = 0;
202:   if (f_max == f_min) *info = 2;
203:   else /* if (complicated condition) */ {
204:     /* Computing MIN */
205:     d__1 = PetscAbsScalar(f_max);
206:     d__2 = PetscAbsScalar(f_min);
207:     if (f_max - f_min > PetscMin(d__1,d__2) * .1) *info = 3;
208:   }
209:   if (*info != 0) return 0;

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

216:   /* Computing MAX */
217:   d__1 = PetscMax(eps[4],eps[5]);
218:   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:   }

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

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

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