1 /* fnoise/snesdnest.F -- translated by f2c (version 20020314). 2 */ 3 #include <petscsys.h> 4 #define FALSE_ 0 5 #define TRUE_ 1 6 7 /* Noise estimation routine, written by Jorge More'. Details are below. */ 8 9 PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *); 10 11 PetscErrorCode SNESNoise_dnest_(PetscInt *nf, double *fval, double *h__, double *fnoise, double *fder2, double *hopt, PetscInt *info, double *eps) 12 { 13 /* Initialized data */ 14 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}; 16 17 /* System generated locals */ 18 PetscInt i__1; 19 double d__1, d__2, d__3, d__4; 20 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; 30 31 /* ********** */ 32 33 /* Subroutine dnest */ 34 35 /* This subroutine estimates the noise in a function */ 36 /* and provides estimates of the optimal difference parameter */ 37 /* for a forward-difference approximation. */ 38 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 */ 42 43 /* f(x-2*h), f(x-h), f(x), f(x+h), f(x+2*h), */ 44 45 /* in the array fval. The use of nf = 7 function evaluations is */ 46 /* recommended. */ 47 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. */ 52 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. */ 58 59 /* The subroutine statement is */ 60 61 /* subroutine dnest(nf,fval,h,hopt,fnoise,info,eps) */ 62 63 /* where */ 64 65 /* nf is a PetscInt variable. */ 66 /* On entry nf is the number of function values. */ 67 /* On exit nf is unchanged. */ 68 69 /* f is a double precision array of dimension nf. */ 70 /* On entry f contains the function values. */ 71 /* On exit f is overwritten. */ 72 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. */ 76 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. */ 81 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. */ 86 87 /* info is a PetscInt variable. */ 88 /* On entry info need not be specified. */ 89 /* On exit info is set as follows: */ 90 91 /* info = 1 Noise has been detected. */ 92 93 /* info = 2 Noise has not been detected; h is too small. */ 94 /* Try 100*h for the next value of h. */ 95 96 /* info = 3 Noise has not been detected; h is too large. */ 97 /* Try h/100 for the next value of h. */ 98 99 /* info = 4 Noise has been detected but the estimate of hopt */ 100 /* is not reliable; h is too small. */ 101 102 /* eps is a double precision work array of dimension nf. */ 103 104 /* MINPACK-2 Project. April 1997. */ 105 /* Argonne National Laboratory. */ 106 /* Jorge J. More'. */ 107 108 /* ********** */ 109 /* Parameter adjustments */ 110 --eps; 111 --fval; 112 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.; 137 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); 147 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_; 163 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 } 170 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); 208 209 /* Determine the noise level. */ 210 /* Computing MIN */ 211 d__1 = PetscMin(eps[4], eps[5]); 212 emin = PetscMin(d__1, eps[6]); 213 214 /* Computing MAX */ 215 d__1 = PetscMax(eps[4], eps[5]); 216 emax = PetscMax(d__1, eps[6]); 217 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 } 229 230 /* Computing MIN */ 231 d__1 = PetscMin(eps[3], eps[4]); 232 emin = PetscMin(d__1, eps[5]); 233 234 /* Computing MAX */ 235 d__1 = PetscMax(eps[3], eps[4]); 236 emax = PetscMax(d__1, eps[5]); 237 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_ */ 270