1 2 /* fnoise/snesdnest.F -- translated by f2c (version 20020314). 3 */ 4 #include <petscsys.h> 5 #define FALSE_ 0 6 #define TRUE_ 1 7 8 /* Noise estimation routine, written by Jorge More'. Details are below. */ 9 10 /* Subroutine */ PetscErrorCode SNESNoise_dnest_(PetscInt *nf, double *fval,double *h__,double *fnoise, double *fder2, double *hopt, PetscInt *info, double *eps) 11 { 12 /* Initialized data */ 13 14 static double const__[15] = { .71,.41,.23,.12,.063,.033,.018,.0089, 15 .0046,.0024,.0012,6.1e-4,3.1e-4,1.6e-4,8e-5 }; 16 17 /* System generated locals */ 18 PetscInt i__1; 19 double d__1, d__2, d__3, d__4; 20 21 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; 31 32 /* ********** */ 33 34 /* Subroutine dnest */ 35 36 /* This subroutine estimates the noise in a function */ 37 /* and provides estimates of the optimal difference parameter */ 38 /* for a forward-difference approximation. */ 39 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 */ 43 44 /* f(x-2*h), f(x-h), f(x), f(x+h), f(x+2*h), */ 45 46 /* in the array fval. The use of nf = 7 function evaluations is */ 47 /* recommended. */ 48 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. */ 53 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 ouput value of info. */ 57 /* In most cases, the subroutine detects noise with the initial */ 58 /* value of h. */ 59 60 /* The subroutine statement is */ 61 62 /* subroutine dnest(nf,fval,h,hopt,fnoise,info,eps) */ 63 64 /* where */ 65 66 /* nf is an int variable. */ 67 /* On entry nf is the number of function values. */ 68 /* On exit nf is unchanged. */ 69 70 /* f is a double precision array of dimension nf. */ 71 /* On entry f contains the function values. */ 72 /* On exit f is overwritten. */ 73 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. */ 77 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. */ 82 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. */ 87 88 /* info is an int variable. */ 89 /* On entry info need not be specified. */ 90 /* On exit info is set as follows: */ 91 92 /* info = 1 Noise has been detected. */ 93 94 /* info = 2 Noise has not been detected; h is too small. */ 95 /* Try 100*h for the next value of h. */ 96 97 /* info = 3 Noise has not been detected; h is too large. */ 98 /* Try h/100 for the next value of h. */ 99 100 /* info = 4 Noise has been detected but the estimate of hopt */ 101 /* is not reliable; h is too small. */ 102 103 /* eps is a double precision work array of dimension nf. */ 104 105 /* MINPACK-2 Project. April 1997. */ 106 /* Argonne National Laboratory. */ 107 /* Jorge J. More'. */ 108 109 /* ********** */ 110 /* Parameter adjustments */ 111 --eps; 112 --fval; 113 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, 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) { 135 *fder2 = est4; 136 } else if (err2 < PetscAbsScalar(est4)) { 137 *fder2 = est3; 138 } else { 139 *fder2 = 0.; 140 } 141 /* Compute the range of function values. */ 142 f_min = fval[1]; 143 f_max = fval[1]; 144 i__1 = *nf; 145 for (i__ = 2; i__ <= i__1; ++i__) { 146 /* Computing MIN */ 147 d__1 = f_min, d__2 = fval[i__]; 148 f_min = PetscMin(d__1,d__2); 149 /* Computing MAX */ 150 d__1 = f_max, 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.) { 163 cancel[j - 1] = TRUE_; 164 } 165 /* Computing MAX */ 166 d__2 = scale, d__3 = (d__1 = fval[i__], PetscAbsScalar(d__1)); 167 scale = PetscMax(d__2,d__3); 168 } 169 /* Compute the estimates for the noise level. */ 170 if (scale == 0.) { 171 stdv = 0.; 172 } else { 173 stdv = 0.; 174 i__1 = *nf - j; 175 for (i__ = 1; i__ <= i__1; ++i__) { 176 /* Computing 2nd power */ 177 d__1 = fval[i__] / scale; 178 stdv += d__1 * d__1; 179 } 180 stdv = scale * PetscSqrtScalar(stdv / (*nf - j)); 181 } 182 eps[j] = const__[j - 1] * stdv; 183 /* Determine differences in sign. */ 184 i__1 = *nf - j - 1; 185 for (i__ = 1; i__ <= i__1; ++i__) { 186 /* Computing MIN */ 187 d__1 = fval[i__], d__2 = fval[i__ + 1]; 188 /* Computing MAX */ 189 d__3 = fval[i__], d__4 = fval[i__ + 1]; 190 if (PetscMin(d__1,d__2) < 0. && PetscMax(d__3,d__4) > 0.) { 191 dsgn[j - 1] = TRUE_; 192 } 193 } 194 } 195 /* First requirement for detection of noise. */ 196 dnoise = dsgn[3]; 197 /* Check for h too small or too large. */ 198 *info = 0; 199 if (f_max == f_min) { 200 *info = 2; 201 } else /* if (complicated condition) */ { 202 /* Computing MIN */ 203 d__1 = PetscAbsScalar(f_max), d__2 = PetscAbsScalar(f_min); 204 if (f_max - f_min > PetscMin(d__1,d__2) * .1) { 205 *info = 3; 206 } 207 } 208 if (*info != 0) { 209 PetscFunctionReturn(0); 210 } 211 /* Determine the noise level. */ 212 /* Computing MIN */ 213 d__1 = PetscMin(eps[4],eps[5]); 214 emin = PetscMin(d__1,eps[6]); 215 /* Computing MAX */ 216 d__1 = PetscMax(eps[4],eps[5]); 217 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(0); 228 } 229 /* Computing MIN */ 230 d__1 = PetscMin(eps[3],eps[4]); 231 emin = PetscMin(d__1,eps[5]); 232 /* Computing MAX */ 233 d__1 = PetscMax(eps[3],eps[4]); 234 emax = PetscMax(d__1,eps[5]); 235 if (emax <= emin * 4 && dnoise) { 236 *fnoise = (eps[3] + eps[4] + eps[5]) / 3; 237 if (*fder2 != 0.) { 238 *info = 1; 239 *hopt = PetscSqrtScalar(*fnoise / PetscAbsScalar(*fder2)) * 1.68; 240 } else { 241 *info = 4; 242 *hopt = *h__ * 10; 243 } 244 PetscFunctionReturn(0); 245 } 246 /* Noise not detected; decide if h is too small or too large. */ 247 if (! cancel[3]) { 248 if (dsgn[3]) { 249 *info = 2; 250 } else { 251 *info = 3; 252 } 253 PetscFunctionReturn(0); 254 } 255 if (! cancel[2]) { 256 if (dsgn[2]) { 257 *info = 2; 258 } else { 259 *info = 3; 260 } 261 PetscFunctionReturn(0); 262 } 263 /* If there is cancelllation on the third and fourth column */ 264 /* then h is too small */ 265 *info = 2; 266 PetscFunctionReturn(0); 267 /* if (cancel .or. dsgn(3)) then */ 268 /* info = 2 */ 269 /* else */ 270 /* info = 3 */ 271 /* end if */ 272 } /* dnest_ */ 273 274