htswanalysis/MACS/lib/gsl/gsl-1.11/min/bracketing.c

   1 /* min/bracketing.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* bracketing.c -- find an initial bracketing interval for a function to minimize */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_min.h>
  26 #include <gsl/gsl_machine.h>
  27
  28 #include "min.h"
  29
  30 int
  31 gsl_min_find_bracket(gsl_function *f,double *x_minimum,double * f_minimum,
  32                      double * x_lower, double * f_lower,
  33                      double * x_upper, double * f_upper,
  34                      size_t eval_max)
  35 {
  36   /* The three following variables must be declared volatile to avoid storage
  37      in extended precision registers available on some architecture. The code
  38      relies on the ability to compare double values. As the values will be
  39      store in regular memory, the extended precision will then be lost and
  40      values that are different in extended precision might have equal
  41      representation in double precision. This behavior might break the
  42      algorithm.
  43    */
  44   volatile double f_left = *f_lower;
  45   volatile double f_right = *f_upper;
  46   volatile double f_center;
  47   double x_left = *x_lower;
  48   double x_right= *x_upper;
  49   double x_center;
  50   const double golden = 0.3819660;      /* golden = (3 - sqrt(5))/2 */
  51   size_t nb_eval = 0;
  52
  53
  54   if (f_right >= f_left)
  55     {
  56       x_center = (x_right - x_left) * golden + x_left;
  57       nb_eval++;
  58       SAFE_FUNC_CALL (f, x_center, &f_center);
  59     }
  60   else
  61     {
  62       x_center = x_right ;
  63       f_center = f_right ;
  64       x_right = (x_center - x_left) / golden + x_left;
  65       nb_eval++;
  66       SAFE_FUNC_CALL (f, x_right, &f_right);
  67     }
  68
  69   do
  70     {
  71       if (f_center < f_left )
  72         {
  73           if (f_center < f_right)
  74             {
  75               *x_lower = x_left;
  76               *x_upper = x_right;
  77               *x_minimum = x_center;
  78               *f_lower = f_left;
  79               *f_upper = f_right;
  80               *f_minimum = f_center;
  81 /*            gsl_ieee_printf_double (&f_left);
  82               printf(" ");
  83               gsl_ieee_printf_double (&f_center);
  84               printf("\n");*/
  85               return GSL_SUCCESS;
  86             }
  87           else if (f_center > f_right)
  88             {
  89               x_left = x_center;
  90               f_left = f_center;
  91               x_center = x_right;
  92               f_center = f_right;
  93               x_right = (x_center - x_left) / golden + x_left;
  94               nb_eval++;
  95               SAFE_FUNC_CALL (f, x_right, &f_right);
  96             }
  97           else /* f_center == f_right */
  98             {
  99               x_right = x_center;
 100               f_right = f_center;
 101               x_center = (x_right - x_left) * golden + x_left;
 102               nb_eval++;
 103               SAFE_FUNC_CALL (f, x_center, &f_center);
 104             }
 105         }
 106       else /* f_center >= f_left */
 107         {
 108           x_right = x_center;
 109           f_right = f_center;
 110           x_center = (x_right - x_left) * golden + x_left;
 111           nb_eval++;
 112           SAFE_FUNC_CALL (f, x_center, &f_center);
 113         }
 114     }
 115   while (nb_eval < eval_max
 116          && (x_right - x_left) > GSL_SQRT_DBL_EPSILON * ( (x_right + x_left) * 0.5 ) + GSL_SQRT_DBL_EPSILON);
 117   *x_lower = x_left;
 118   *x_upper = x_right;
 119   *x_minimum = x_center;
 120   *f_lower = f_left;
 121   *f_upper = f_right;
 122   *f_minimum = f_center;
 123   return GSL_FAILURE;
 124 }
 125