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_ */