/* multimin/conjugate_fr.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. *//* conjugate_fr.c -- Conjugate gradient Fletcher-Reeve algorithm *//* Modified by Brian Gough to use single iteration structure */#include <config.h>#include <gsl/gsl_multimin.h>#include <gsl/gsl_blas.h>#include "directional_minimize.c"typedef struct{  int iter;  double step;  double max_step;  double tol;  gsl_vector *x1;  gsl_vector *dx1;  gsl_vector *x2;  double pnorm;  gsl_vector *p;  double g0norm;  gsl_vector *g0;}conjugate_fr_state_t;static intconjugate_fr_alloc (void *vstate, size_t n){  conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;  state->x1 = gsl_vector_calloc (n);  if (state->x1 == 0)    {      GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);    }  state->dx1 = gsl_vector_calloc (n);  if (state->dx1 == 0)    {      gsl_vector_free (state->x1);      GSL_ERROR ("failed to allocate space for dx1", GSL_ENOMEM);    }  state->x2 = gsl_vector_calloc (n);  if (state->x2 == 0)    {      gsl_vector_free (state->dx1);      gsl_vector_free (state->x1);      GSL_ERROR ("failed to allocate space for x2", GSL_ENOMEM);    }  state->p = gsl_vector_calloc (n);  if (state->p == 0)    {      gsl_vector_free (state->x2);      gsl_vector_free (state->dx1);      gsl_vector_free (state->x1);      GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);    }  state->g0 = gsl_vector_calloc (n);  if (state->g0 == 0)    {      gsl_vector_free (state->p);      gsl_vector_free (state->x2);      gsl_vector_free (state->dx1);      gsl_vector_free (state->x1);      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);    }  return GSL_SUCCESS;}static intconjugate_fr_set (void *vstate, gsl_multimin_function_fdf * fdf,                  const gsl_vector * x, double *f, gsl_vector * gradient,                  double step_size, double tol){  conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;  state->iter = 0;  state->step = step_size;  state->max_step = step_size;  state->tol  = tol;  GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);  /* Use the gradient as the initial direction */  gsl_vector_memcpy (state->p, gradient);  gsl_vector_memcpy (state->g0, gradient);  {    double gnorm = gsl_blas_dnrm2 (gradient);    state->pnorm = gnorm;    state->g0norm = gnorm;  }  return GSL_SUCCESS;}static voidconjugate_fr_free (void *vstate){  conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;  gsl_vector_free (state->g0);  gsl_vector_free (state->p);  gsl_vector_free (state->x2);  gsl_vector_free (state->dx1);  gsl_vector_free (state->x1);}static intconjugate_fr_restart (void *vstate){  conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;  state->iter = 0;  return GSL_SUCCESS;}static intconjugate_fr_iterate (void *vstate, gsl_multimin_function_fdf * fdf,                      gsl_vector * x, double *f,                      gsl_vector * gradient, gsl_vector * dx){  conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;  gsl_vector *x1 = state->x1;  gsl_vector *dx1 = state->dx1;  gsl_vector *x2 = state->x2;  gsl_vector *p = state->p;  gsl_vector *g0 = state->g0;  double pnorm = state->pnorm;  double g0norm = state->g0norm;  double fa = *f, fb, fc;  double dir;  double stepa = 0.0, stepb, stepc = state->step, tol = state->tol;  double g1norm;  double pg;  if (pnorm == 0.0 || g0norm == 0.0)    {      gsl_vector_set_zero (dx);      return GSL_ENOPROG;    }    /* Determine which direction is downhill, +p or -p */  gsl_blas_ddot (p, gradient, &pg);  dir = (pg >= 0.0) ? +1.0 : -1.0;  /* Compute new trial point at x_c= x - step * p, where p is the     current direction */  take_step (x, p, stepc, dir / pnorm, x1, dx);  /* Evaluate function and gradient at new point xc */  fc = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);  if (fc < fa)    {      /* Success, reduced the function value */      state->step = stepc * 2.0;      *f = fc;      gsl_vector_memcpy (x, x1);      GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);      return GSL_SUCCESS;    }#ifdef DEBUG  printf ("got stepc = %g fc = %g\n", stepc, fc);#endif  /* Do a line minimisation in the region (xa,fa) (xc,fc) to find an     intermediate (xb,fb) satisifying fa > fb < fc.  Choose an initial     xb based on parabolic interpolation */  intermediate_point (fdf, x, p, dir / pnorm, pg,                      stepa, stepc, fa, fc, x1, dx1, gradient, &stepb, &fb);  if (stepb == 0.0)    {      return GSL_ENOPROG;    }  minimize (fdf, x, p, dir / pnorm,            stepa, stepb, stepc, fa, fb, fc, tol,            x1, dx1, x2, dx, gradient, &(state->step), f, &g1norm);  gsl_vector_memcpy (x, x2);  /* Choose a new conjugate direction for the next step */  state->iter = (state->iter + 1) % x->size;  if (state->iter == 0)    {      gsl_vector_memcpy (p, gradient);      state->pnorm = g1norm;    }  else    {      /* p' = g1 - beta * p */      double beta = -pow (g1norm / g0norm, 2.0);      gsl_blas_dscal (-beta, p);      gsl_blas_daxpy (1.0, gradient, p);      state->pnorm = gsl_blas_dnrm2 (p);    }  state->g0norm = g1norm;  gsl_vector_memcpy (g0, gradient);#ifdef DEBUG  printf ("updated conjugate directions\n");  printf ("p: ");  gsl_vector_fprintf (stdout, p, "%g");  printf ("g: ");  gsl_vector_fprintf (stdout, gradient, "%g");#endif  return GSL_SUCCESS;}static const gsl_multimin_fdfminimizer_type conjugate_fr_type = {  "conjugate_fr",               /* name */  sizeof (conjugate_fr_state_t),  &conjugate_fr_alloc,  &conjugate_fr_set,  &conjugate_fr_iterate,  &conjugate_fr_restart,  &conjugate_fr_free};const gsl_multimin_fdfminimizer_type  * gsl_multimin_fdfminimizer_conjugate_fr = &conjugate_fr_type;