1 files changed, 171 insertions, 0 deletions
diff --git a/gsl-1.9/multimin/steepest_descent.c b/gsl-1.9/multimin/steepest_descent.c
new file mode 100644
index 0000000..ebfd29b
--- /dev/null
+++ b/gsl-1.9/multimin/steepest_descent.c
@@ -0,0 +1,171 @@
+/* multimin/steepest_descent.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.
+ */
+
+/* steepest_descent.c -- the steepest descent algorithm */
+
+/* Modified by Brian Gough to use single iteration structure */
+
+#include <config.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_blas_types.h>
+#include <gsl/gsl_blas.h>
+
+typedef struct
+{
+  double step;
+  double max_step;
+  double tol;
+  gsl_vector *x1;
+  gsl_vector *g1;
+}
+steepest_descent_state_t;
+
+static int
+steepest_descent_alloc (void *vstate, size_t n)
+{
+  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
+
+  state->x1 = gsl_vector_alloc (n);
+
+  if (state->x1 == NULL)
+    {
+      GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
+    }
+
+  state->g1 = gsl_vector_alloc (n);
+
+  if (state->g1 == NULL)
+    {
+      gsl_vector_free (state->x1);
+      GSL_ERROR ("failed to allocate space for g1", GSL_ENOMEM);
+    }
+
+  return GSL_SUCCESS;
+}
+
+static int
+steepest_descent_set (void *vstate, gsl_multimin_function_fdf * fdf,
+                      const gsl_vector * x, double *f,
+                      gsl_vector * gradient, double step_size, double tol)
+{
+  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
+
+  GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);
+
+  state->step = step_size;
+  state->max_step = step_size;
+  state->tol = tol;
+
+  return GSL_SUCCESS;
+}
+
+
+static void
+steepest_descent_free (void *vstate)
+{
+  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
+
+  gsl_vector_free (state->x1);
+  gsl_vector_free (state->g1);
+}
+
+static int
+steepest_descent_restart (void *vstate)
+{
+  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
+
+  state->step = state->max_step;
+
+  return GSL_SUCCESS;
+}
+
+static int
+steepest_descent_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
+                          gsl_vector * x, double *f,
+                          gsl_vector * gradient, gsl_vector * dx)
+{
+  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;
+
+  gsl_vector *x1 = state->x1;
+  gsl_vector *g1 = state->g1;
+
+  double f0 = *f;
+  double f1;
+  double step = state->step, tol = state->tol;
+
+  int failed = 0;
+
+  /* compute new trial point at x1= x - step * dir, where dir is the
+     normalized gradient */
+
+  double gnorm = gsl_blas_dnrm2 (gradient);
+
+  if (gnorm == 0.0)
+    {
+      gsl_vector_set_zero (dx);
+      return GSL_ENOPROG;
+    }
+
+trial:
+  gsl_vector_set_zero (dx);
+  gsl_blas_daxpy (-step / gnorm, gradient, dx);
+
+  gsl_vector_memcpy (x1, x);
+  gsl_blas_daxpy (1.0, dx, x1);
+
+  /* evaluate function and gradient at new point x1 */
+
+  GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x1, &f1, g1);
+
+  if (f1 > f0)
+    {
+      /* downhill step failed, reduce step-size and try again */
+
+      failed = 1;
+      step *= tol;
+      goto trial;
+    }
+
+  if (failed)
+    step *= tol;
+  else
+    step *= 2.0;
+
+  state->step = step;
+
+  gsl_vector_memcpy (x, x1);
+  gsl_vector_memcpy (gradient, g1);
+
+  *f = f1;
+
+  return GSL_SUCCESS;
+}
+
+static const gsl_multimin_fdfminimizer_type steepest_descent_type =
+  { "steepest_descent",         /* name */
+  sizeof (steepest_descent_state_t),
+  &steepest_descent_alloc,
+  &steepest_descent_set,
+  &steepest_descent_iterate,
+  &steepest_descent_restart,
+  &steepest_descent_free
+};
+
+const gsl_multimin_fdfminimizer_type
+  * gsl_multimin_fdfminimizer_steepest_descent = &steepest_descent_type;