typedef struct

{

  gsl_function_fdf fdf_linear;

  gsl_multimin_function_fdf *fdf;

  /* fixed values */

  const gsl_vector *x;

  const gsl_vector *g;

  const gsl_vector *p;

  /* cached values, for x(alpha) = x + alpha * p */

  double f_alpha;

  double df_alpha;

  gsl_vector *x_alpha;

  gsl_vector *g_alpha;

  /* cache "keys" */

  double f_cache_key;

  double df_cache_key;

  double x_cache_key;

  double g_cache_key;

}

wrapper_t;

static void

moveto (double alpha, wrapper_t * w)

{

  if (alpha == w->x_cache_key)  /* using previously cached position */

    {

      return;

    }

  /* set x_alpha = x + alpha * p */

  gsl_vector_memcpy (w->x_alpha, w->x);

  gsl_blas_daxpy (alpha, w->p, w->x_alpha);

  w->x_cache_key = alpha;

}

static double

slope (wrapper_t * w)           /* compute gradient . direction */

{

  double df;

  gsl_blas_ddot (w->g_alpha, w->p, &df);

  return df;

}

static double

wrap_f (double alpha, void *params)

{

  wrapper_t *w = (wrapper_t *) params;

  if (alpha == w->f_cache_key)  /* using previously cached f(alpha) */

    {

      return w->f_alpha;

    }

  moveto (alpha, w);

  w->f_alpha = GSL_MULTIMIN_FN_EVAL_F (w->fdf, w->x_alpha);

  w->f_cache_key = alpha;

  return w->f_alpha;

}

static double

wrap_df (double alpha, void *params)

{

  wrapper_t *w = (wrapper_t *) params;

  if (alpha == w->df_cache_key) /* using previously cached df(alpha) */

    {

      return w->df_alpha;

    }

  moveto (alpha, w);

  if (alpha != w->g_cache_key) 

    {

      GSL_MULTIMIN_FN_EVAL_DF (w->fdf, w->x_alpha, w->g_alpha);

      w->g_cache_key = alpha;

    }

  w->df_alpha = slope (w);

  w->df_cache_key = alpha;

  return w->df_alpha;

}

static void

wrap_fdf (double alpha, void *params, double *f, double *df)

{

  wrapper_t *w = (wrapper_t *) params;

  /* Check for previously cached values */

  if (alpha == w->f_cache_key && alpha == w->df_cache_key)

    {

      *f = w->f_alpha;

      *df = w->df_alpha;

      return;

    }

  if (alpha == w->f_cache_key || alpha == w->df_cache_key)

    {

      *f = wrap_f (alpha, params);

      *df = wrap_df (alpha, params);

      return;

    }

  moveto (alpha, w);

  GSL_MULTIMIN_FN_EVAL_F_DF (w->fdf, w->x_alpha, &w->f_alpha, w->g_alpha);

  w->f_cache_key = alpha;

  w->g_cache_key = alpha;

  w->df_alpha = slope (w);

  w->df_cache_key = alpha;

  *f = w->f_alpha;

  *df = w->df_alpha;

}

static void

prepare_wrapper (wrapper_t * w, gsl_multimin_function_fdf * fdf,

                 const gsl_vector * x, double f, const gsl_vector *g, 

                 const gsl_vector * p,

                 gsl_vector * x_alpha, gsl_vector *g_alpha)

{

  w->fdf_linear.f = &wrap_f;

  w->fdf_linear.df = &wrap_df;

  w->fdf_linear.fdf = &wrap_fdf;

  w->fdf_linear.params = (void *)w;  /* pointer to "self" */

  w->fdf = fdf;

  w->x = x;

  w->g = g;

  w->p = p;

  w->x_alpha = x_alpha;

  w->g_alpha = g_alpha;

  gsl_vector_memcpy(w->x_alpha, w->x); 

  w->x_cache_key = 0.0;

  w->f_alpha = f;

  w->f_cache_key = 0.0;

  gsl_vector_memcpy(w->g_alpha, w->g);

  w->g_cache_key = 0.0;

  w->df_alpha = slope(w);

  w->df_cache_key = 0.0;

}

static void

update_position (wrapper_t * w, double alpha, gsl_vector *x, double *f, gsl_vector *g)

{

  /* ensure that everything is fully cached */

  { double f_alpha, df_alpha; wrap_fdf (alpha, w, &f_alpha, &df_alpha); } ;

  *f = w->f_alpha;

  gsl_vector_memcpy(x, w->x_alpha);

  gsl_vector_memcpy(g, w->g_alpha);

}  

static void

change_direction (wrapper_t * w)

{

  /* Convert the cache values from the end of the current minimisation

     to those needed for the start of the next minimisation, alpha=0 */

  /* The new x_alpha for alpha=0 is the current position */

  gsl_vector_memcpy (w->x_alpha, w->x);

  w->x_cache_key = 0.0;

  /* The function value does not change */

  w->f_cache_key = 0.0;

  /* The new g_alpha for alpha=0 is the current gradient at the endpoint */

  gsl_vector_memcpy (w->g_alpha, w->g);

  w->g_cache_key = 0.0;

  /* Calculate the slope along the new direction vector, p */

  w->df_alpha = slope (w);

  w->df_cache_key = 0.0;

}