Actual source code: snesdnest.c

      1: #define PETSCSNES_DLL

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

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

     11: /* Subroutine */ PetscErrorCode 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,
     16:             .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;

     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 an int 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 an int 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__ *
    124:              2);
    125:     est3 = (fval[mh + 3] - fval[mh] * 2 + fval[mh - 3]) / (*h__ * 3) / (*h__ *
    126:              3);
    127:     est4 = (est1 + est2 + est3) / 3;
    128: /* Computing MAX */
    129: /* Computing PETSCMAX */
    130:     d__3 = PetscMax(est1,est2);
    131: /* Computing MIN */
    132:     d__4 = PetscMin(est1,est2);
    133:     d__1 = PetscMax(d__3,est3) - est4, d__2 = est4 - PetscMin(d__4,est3);
    134:     err2 = PetscMax(d__1,d__2);
    135: /*      write (2,123) est1, est2, est3 */
    136: /* 123  format ('Second derivative estimates', 3d12.2) */
    137:     if (err2 <= PetscAbsScalar(est4) * .1) {
    138:         *fder2 = est4;
    139:     } else if (err2 < PetscAbsScalar(est4)) {
    140:         *fder2 = est3;
    141:     } else {
    142:         *fder2 = 0.;
    143:     }
    144: /*     Compute the range of function values. */
    145:     f_min = fval[1];
    146:     f_max = fval[1];
    147:     i__1 = *nf;
    148:     for (i__ = 2; i__ <= i__1; ++i__) {
    149: /* Computing MIN */
    150:         d__1 = f_min, d__2 = fval[i__];
    151:         f_min = PetscMin(d__1,d__2);
    152: /* Computing MAX */
    153:         d__1 = f_max, d__2 = fval[i__];
    154:         f_max = PetscMax(d__1,d__2);
    155:     }
    156: /*     Construct the difference table. */
    157:     dnoise = FALSE_;
    158:     for (j = 1; j <= 6; ++j) {
    159:         dsgn[j - 1] = FALSE_;
    160:         cancel[j - 1] = FALSE_;
    161:         scale = 0.;
    162:         i__1 = *nf - j;
    163:         for (i__ = 1; i__ <= i__1; ++i__) {
    164:             fval[i__] = fval[i__ + 1] - fval[i__];
    165:             if (fval[i__] == 0.) {
    166:                 cancel[j - 1] = TRUE_;
    167:             }
    168: /* Computing MAX */
    169:             d__2 = scale, d__3 = (d__1 = fval[i__], PetscAbsScalar(d__1));
    170:             scale = PetscMax(d__2,d__3);
    171:         }
    172: /*        Compute the estimates for the noise level. */
    173:         if (scale == 0.) {
    174:             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__], d__2 = fval[i__ + 1];
    191: /* Computing MAX */
    192:             d__3 = fval[i__], d__4 = fval[i__ + 1];
    193:             if (PetscMin(d__1,d__2) < 0. && PetscMax(d__3,d__4) > 0.) {
    194:                 dsgn[j - 1] = TRUE_;
    195:             }
    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) {
    203:         *info = 2;
    204:     } else /* if(complicated condition) */ {
    205: /* Computing MIN */
    206:         d__1 = PetscAbsScalar(f_max), d__2 = PetscAbsScalar(f_min);
    207:         if (f_max - f_min > PetscMin(d__1,d__2) * .1) {
    208:             *info = 3;
    209:         }
    210:     }
    211:     if (*info != 0) {
    212:         return(0);
    213:     }
    214: /*     Determine the noise level. */
    215: /* Computing MIN */
    216:     d__1 = PetscMin(eps[4],eps[5]);
    217:     emin = PetscMin(d__1,eps[6]);
    218: /* Computing MAX */
    219:     d__1 = PetscMax(eps[4],eps[5]);
    220:     emax = PetscMax(d__1,eps[6]);
    221:     if (emax <= emin * 4 && dnoise) {
    222:         *fnoise = (eps[4] + eps[5] + eps[6]) / 3;
    223:         if (*fder2 != 0.) {
    224:             *info = 1;
    225:             *hopt = PetscSqrtScalar(*fnoise / PetscAbsScalar(*fder2)) * 1.68;
    226:         } else {
    227:             *info = 4;
    228:             *hopt = *h__ * 10;
    229:         }
    230:         return(0);
    231:     }
    232: /* Computing MIN */
    233:     d__1 = PetscMin(eps[3],eps[4]);
    234:     emin = PetscMin(d__1,eps[5]);
    235: /* Computing MAX */
    236:     d__1 = PetscMax(eps[3],eps[4]);
    237:     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:         return(0);
    248:     }
    249: /*     Noise not detected; decide if h is too small or too large. */
    250:     if (! cancel[3]) {
    251:         if (dsgn[3]) {
    252:             *info = 2;
    253:         } else {
    254:             *info = 3;
    255:         }
    256:         return(0);
    257:     }
    258:     if (! cancel[2]) {
    259:         if (dsgn[2]) {
    260:             *info = 2;
    261:         } else {
    262:             *info = 3;
    263:         }
    264:         return(0);
    265:     }
    266: /*     If there is cancelllation on the third and fourth column */
    267: /*     then h is too small */
    268:     *info = 2;
    269:     return(0);
    270: /*      if (cancel .or. dsgn(3)) then */
    271: /*         info = 2 */
    272: /*      else */
    273: /*         info = 3 */
    274: /*      end if */
    275: } /* dnest_ */