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

   1 /* multimin/steepest_descent.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 /* steepest_descent.c -- the steepest descent algorithm */
  21
  22 /* Modified by Brian Gough to use single iteration structure */
  23
  24 #include <config.h>
  25 #include <gsl/gsl_multimin.h>
  26 #include <gsl/gsl_blas_types.h>
  27 #include <gsl/gsl_blas.h>
  28
  29 typedef struct
  30 {
  31   double step;
  32   double max_step;
  33   double tol;
  34   gsl_vector *x1;
  35   gsl_vector *g1;
  36 }
  37 steepest_descent_state_t;
  38
  39 static int
  40 steepest_descent_alloc (void *vstate, size_t n)
  41 {
  42   steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
  43
  44   state->x1 = gsl_vector_alloc (n);
  45
  46   if (state->x1 == NULL)
  47     {
  48       GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
  49     }
  50
  51   state->g1 = gsl_vector_alloc (n);
  52
  53   if (state->g1 == NULL)
  54     {
  55       gsl_vector_free (state->x1);
  56       GSL_ERROR ("failed to allocate space for g1", GSL_ENOMEM);
  57     }
  58
  59   return GSL_SUCCESS;
  60 }
  61
  62 static int
  63 steepest_descent_set (void *vstate, gsl_multimin_function_fdf * fdf,
  64                       const gsl_vector * x, double *f,
  65                       gsl_vector * gradient, double step_size, double tol)
  66 {
  67   steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
  68
  69   GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);
  70
  71   state->step = step_size;
  72   state->max_step = step_size;
  73   state->tol = tol;
  74
  75   return GSL_SUCCESS;
  76 }
  77
  78
  79 static void
  80 steepest_descent_free (void *vstate)
  81 {
  82   steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
  83
  84   gsl_vector_free (state->x1);
  85   gsl_vector_free (state->g1);
  86 }
  87
  88 static int
  89 steepest_descent_restart (void *vstate)
  90 {
  91   steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
  92
  93   state->step = state->max_step;
  94
  95   return GSL_SUCCESS;
  96 }
  97
  98 static int
  99 steepest_descent_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
 100                           gsl_vector * x, double *f,
 101                           gsl_vector * gradient, gsl_vector * dx)
 102 {
 103   steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
 104
 105   gsl_vector *x1 = state->x1;
 106   gsl_vector *g1 = state->g1;
 107
 108   double f0 = *f;
 109   double f1;
 110   double step = state->step, tol = state->tol;
 111
 112   int failed = 0;
 113
 114   /* compute new trial point at x1= x - step * dir, where dir is the
 115      normalized gradient */
 116
 117   double gnorm = gsl_blas_dnrm2 (gradient);
 118
 119   if (gnorm == 0.0)
 120     {
 121       gsl_vector_set_zero (dx);
 122       return GSL_ENOPROG;
 123     }
 124
 125 trial:
 126   gsl_vector_set_zero (dx);
 127   gsl_blas_daxpy (-step / gnorm, gradient, dx);
 128
 129   gsl_vector_memcpy (x1, x);
 130   gsl_blas_daxpy (1.0, dx, x1);
 131
 132   /* evaluate function and gradient at new point x1 */
 133
 134   GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x1, &f1, g1);
 135
 136   if (f1 > f0)
 137     {
 138       /* downhill step failed, reduce step-size and try again */
 139
 140       failed = 1;
 141       step *= tol;
 142       goto trial;
 143     }
 144
 145   if (failed)
 146     step *= tol;
 147   else
 148     step *= 2.0;
 149
 150   state->step = step;
 151
 152   gsl_vector_memcpy (x, x1);
 153   gsl_vector_memcpy (gradient, g1);
 154
 155   *f = f1;
 156
 157   return GSL_SUCCESS;
 158 }
 159
 160 static const gsl_multimin_fdfminimizer_type steepest_descent_type =
 161   { "steepest_descent",         /* name */
 162   sizeof (steepest_descent_state_t),
 163   &steepest_descent_alloc,
 164   &steepest_descent_set,
 165   &steepest_descent_iterate,
 166   &steepest_descent_restart,
 167   &steepest_descent_free
 168 };
 169
 170 const gsl_multimin_fdfminimizer_type
 171   * gsl_multimin_fdfminimizer_steepest_descent = &steepest_descent_type;