htswanalysis/MACS/lib/gsl/gsl-1.11/multimin/diff.c

   1 /* multimin/diff.c
   2  *
   3  * Copyright (C) 2000 David Morrison
   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 #include <config.h>
  21 #include <gsl/gsl_multimin.h>
  22
  23 int
  24 gsl_multimin_diff (const gsl_multimin_function * f,
  25                    const gsl_vector * x, gsl_vector * g)
  26 {
  27   size_t i, n = f->n;
  28
  29   double h = GSL_SQRT_DBL_EPSILON;
  30
  31
  32   gsl_vector * x1 = gsl_vector_alloc (n);  /* FIXME: pass as argument */
  33
  34   gsl_vector_memcpy (x1, x);
  35
  36   for (i = 0; i < n; i++)
  37     {
  38       double fl, fh;
  39
  40       double xi = gsl_vector_get (x, i);
  41       double dx = fabs(xi) * h;
  42
  43       if (dx == 0.0) dx = h;
  44
  45       gsl_vector_set (x1, i, xi + dx);
  46       fh = GSL_MULTIMIN_FN_EVAL(f, x1);
  47
  48       gsl_vector_set (x1, i, xi - dx);
  49       fl = GSL_MULTIMIN_FN_EVAL(f, x1);
  50
  51       gsl_vector_set (x1, i, xi);
  52       gsl_vector_set (g, i, (fh - fl) / (2.0 * dx));
  53     }
  54
  55   gsl_vector_free (x1);
  56
  57   return GSL_SUCCESS;
  58 }