diff options
Diffstat (limited to 'gsl-1.9/min/bracketing.c')
| -rw-r--r-- | gsl-1.9/min/bracketing.c | 125 |
1 files changed, 125 insertions, 0 deletions
diff --git a/gsl-1.9/min/bracketing.c b/gsl-1.9/min/bracketing.c new file mode 100644 index 0000000..f3ceaf4 --- /dev/null +++ b/gsl-1.9/min/bracketing.c @@ -0,0 +1,125 @@ +/* min/bracketing.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* bracketing.c -- find an initial bracketing interval for a function to minimize */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_min.h> +#include <gsl/gsl_machine.h> + +#include "min.h" + +int +gsl_min_find_bracket(gsl_function *f,double *x_minimum,double * f_minimum, + double * x_lower, double * f_lower, + double * x_upper, double * f_upper, + size_t eval_max) +{ + /* The three following variables must be declared volatile to avoid storage + in extended precision registers available on some architecture. The code + relies on the ability to compare double values. As the values will be + store in regular memory, the extended precision will then be lost and + values that are different in extended precision might have equal + representation in double precision. This behavior might break the + algorithm. + */ + volatile double f_left = *f_lower; + volatile double f_right = *f_upper; + volatile double f_center; + double x_left = *x_lower; + double x_right= *x_upper; + double x_center; + const double golden = 0.3819660; /* golden = (3 - sqrt(5))/2 */ + size_t nb_eval = 0; + + + if (f_right >= f_left) + { + x_center = (x_right - x_left) * golden + x_left; + nb_eval++; + SAFE_FUNC_CALL (f, x_center, &f_center); + } + else + { + x_center = x_right ; + f_center = f_right ; + x_right = (x_center - x_left) / golden + x_left; + nb_eval++; + SAFE_FUNC_CALL (f, x_right, &f_right); + } + + do + { + if (f_center < f_left ) + { + if (f_center < f_right) + { + *x_lower = x_left; + *x_upper = x_right; + *x_minimum = x_center; + *f_lower = f_left; + *f_upper = f_right; + *f_minimum = f_center; +/* gsl_ieee_printf_double (&f_left); + printf(" "); + gsl_ieee_printf_double (&f_center); + printf("\n");*/ + return GSL_SUCCESS; + } + else if (f_center > f_right) + { + x_left = x_center; + f_left = f_center; + x_center = x_right; + f_center = f_right; + x_right = (x_center - x_left) / golden + x_left; + nb_eval++; + SAFE_FUNC_CALL (f, x_right, &f_right); + } + else /* f_center == f_right */ + { + x_right = x_center; + f_right = f_center; + x_center = (x_right - x_left) * golden + x_left; + nb_eval++; + SAFE_FUNC_CALL (f, x_center, &f_center); + } + } + else /* f_center >= f_left */ + { + x_right = x_center; + f_right = f_center; + x_center = (x_right - x_left) * golden + x_left; + nb_eval++; + SAFE_FUNC_CALL (f, x_center, &f_center); + } + } + while (nb_eval < eval_max + && (x_right - x_left) > GSL_SQRT_DBL_EPSILON * ( (x_right + x_left) * 0.5 ) + GSL_SQRT_DBL_EPSILON); + *x_lower = x_left; + *x_upper = x_right; + *x_minimum = x_center; + *f_lower = f_left; + *f_upper = f_right; + *f_minimum = f_center; + return GSL_FAILURE; +} + |