1 /* fnoise/snesdnest.F -- translated by f2c (version 20020314). 2 */ 3 #include "petsc.h" 4 #define FALSE_ 0 5 #define TRUE_ 1 6 7 /* "$Id: snesdnest.F,v 1.6 1998/03/25 17:15:03 balay Exp bsmith $"; */ 8 9 /* Noise estimation routine, written by Jorge More'. Details are below. */ 10 11 /* Subroutine */ int dnest_(int *nf, double *fval,double *h__,double *fnoise, double *fder2, double *hopt, int *info, double *eps) 12 { 13 /* Initialized data */ 14 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 }; 17 18 /* System generated locals */ 19 int i__1; 20 double d__1, d__2, d__3, d__4; 21 22 23 /* Local variables */ 24 static double emin, emax; 25 static int dsgn[6]; 26 static double fmax, fmin, stdv; 27 static int i__, j; 28 static double scale; 29 static int mh; 30 static int cancel[6], dnoise; 31 static double err2, est1, est2, est3, est4; 32 33 /* ********** */ 34 35 /* Subroutine dnest */ 36 37 /* This subroutine estimates the noise in a function */ 38 /* and provides estimates of the optimal difference parameter */ 39 /* for a forward-difference approximation. */ 40 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 */ 44 45 /* f(x-2*h), f(x-h), f(x), f(x+h), f(x+2*h), */ 46 47 /* in the array fval. The use of nf = 7 function evaluations is */ 48 /* recommended. */ 49 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. */ 54 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. */ 60 61 /* The subroutine statement is */ 62 63 /* subroutine dnest(nf,fval,h,hopt,fnoise,info,eps) */ 64 65 /* where */ 66 67 /* nf is an int variable. */ 68 /* On entry nf is the number of function values. */ 69 /* On exit nf is unchanged. */ 70 71 /* f is a double precision array of dimension nf. */ 72 /* On entry f contains the function values. */ 73 /* On exit f is overwritten. */ 74 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. */ 78 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. */ 83 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. */ 88 89 /* info is an int variable. */ 90 /* On entry info need not be specified. */ 91 /* On exit info is set as follows: */ 92 93 /* info = 1 Noise has been detected. */ 94 95 /* info = 2 Noise has not been detected; h is too small. */ 96 /* Try 100*h for the next value of h. */ 97 98 /* info = 3 Noise has not been detected; h is too large. */ 99 /* Try h/100 for the next value of h. */ 100 101 /* info = 4 Noise has been detected but the estimate of hopt */ 102 /* is not reliable; h is too small. */ 103 104 /* eps is a double precision work array of dimension nf. */ 105 106 /* MINPACK-2 Project. April 1997. */ 107 /* Argonne National Laboratory. */ 108 /* Jorge J. More'. */ 109 110 /* ********** */ 111 /* Parameter adjustments */ 112 --eps; 113 --fval; 114 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 fmin = fval[1]; 146 fmax = fval[1]; 147 i__1 = *nf; 148 for (i__ = 2; i__ <= i__1; ++i__) { 149 /* Computing MIN */ 150 d__1 = fmin, d__2 = fval[i__]; 151 fmin = PetscMin(d__1,d__2); 152 /* Computing MAX */ 153 d__1 = fmax, d__2 = fval[i__]; 154 fmax = 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 (fmax == fmin) { 203 *info = 2; 204 } else /* if(complicated condition) */ { 205 /* Computing MIN */ 206 d__1 = PetscAbsScalar(fmax), d__2 = PetscAbsScalar(fmin); 207 if (fmax - fmin > 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_ */ 276 277