|
Packit |
67cb25 |
/* multimin/conjugate_fr.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* conjugate_fr.c -- Conjugate gradient Fletcher-Reeve algorithm */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Modified by Brian Gough to use single iteration structure */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_multimin.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "directional_minimize.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int iter;
|
|
Packit |
67cb25 |
double step;
|
|
Packit |
67cb25 |
double max_step;
|
|
Packit |
67cb25 |
double tol;
|
|
Packit |
67cb25 |
gsl_vector *x1;
|
|
Packit |
67cb25 |
gsl_vector *dx1;
|
|
Packit |
67cb25 |
gsl_vector *x2;
|
|
Packit |
67cb25 |
double pnorm;
|
|
Packit |
67cb25 |
gsl_vector *p;
|
|
Packit |
67cb25 |
double g0norm;
|
|
Packit |
67cb25 |
gsl_vector *g0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
conjugate_fr_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
conjugate_fr_alloc (void *vstate, size_t n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->x1 = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->x1 == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->dx1 = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->dx1 == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_free (state->x1);
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for dx1", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->x2 = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->x2 == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_free (state->dx1);
|
|
Packit |
67cb25 |
gsl_vector_free (state->x1);
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for x2", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->p = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->p == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_free (state->x2);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dx1);
|
|
Packit |
67cb25 |
gsl_vector_free (state->x1);
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->g0 = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->g0 == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_free (state->p);
|
|
Packit |
67cb25 |
gsl_vector_free (state->x2);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dx1);
|
|
Packit |
67cb25 |
gsl_vector_free (state->x1);
|
|
Packit |
67cb25 |
GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
conjugate_fr_set (void *vstate, gsl_multimin_function_fdf * fdf,
|
|
Packit |
67cb25 |
const gsl_vector * x, double *f, gsl_vector * gradient,
|
|
Packit |
67cb25 |
double step_size, double tol)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->iter = 0;
|
|
Packit |
67cb25 |
state->step = step_size;
|
|
Packit |
67cb25 |
state->max_step = step_size;
|
|
Packit |
67cb25 |
state->tol = tol;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Use the gradient as the initial direction */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_memcpy (state->p, gradient);
|
|
Packit |
67cb25 |
gsl_vector_memcpy (state->g0, gradient);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double gnorm = gsl_blas_dnrm2 (gradient);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->pnorm = gnorm;
|
|
Packit |
67cb25 |
state->g0norm = gnorm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
conjugate_fr_free (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free (state->g0);
|
|
Packit |
67cb25 |
gsl_vector_free (state->p);
|
|
Packit |
67cb25 |
gsl_vector_free (state->x2);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dx1);
|
|
Packit |
67cb25 |
gsl_vector_free (state->x1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
conjugate_fr_restart (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->iter = 0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
conjugate_fr_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
|
|
Packit |
67cb25 |
gsl_vector * x, double *f,
|
|
Packit |
67cb25 |
gsl_vector * gradient, gsl_vector * dx)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
conjugate_fr_state_t *state = (conjugate_fr_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector *x1 = state->x1;
|
|
Packit |
67cb25 |
gsl_vector *dx1 = state->dx1;
|
|
Packit |
67cb25 |
gsl_vector *x2 = state->x2;
|
|
Packit |
67cb25 |
gsl_vector *p = state->p;
|
|
Packit |
67cb25 |
gsl_vector *g0 = state->g0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double pnorm = state->pnorm;
|
|
Packit |
67cb25 |
double g0norm = state->g0norm;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double fa = *f, fb, fc;
|
|
Packit |
67cb25 |
double dir;
|
|
Packit |
67cb25 |
double stepa = 0.0, stepb, stepc = state->step, tol = state->tol;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double g1norm;
|
|
Packit |
67cb25 |
double pg;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (pnorm == 0.0 || g0norm == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_set_zero (dx);
|
|
Packit |
67cb25 |
return GSL_ENOPROG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Determine which direction is downhill, +p or -p */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_ddot (p, gradient, &pg;;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
dir = (pg >= 0.0) ? +1.0 : -1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute new trial point at x_c= x - step * p, where p is the
|
|
Packit |
67cb25 |
current direction */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
take_step (x, p, stepc, dir / pnorm, x1, dx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate function and gradient at new point xc */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fc = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fc < fa)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Success, reduced the function value */
|
|
Packit |
67cb25 |
state->step = stepc * 2.0;
|
|
Packit |
67cb25 |
*f = fc;
|
|
Packit |
67cb25 |
gsl_vector_memcpy (x, x1);
|
|
Packit |
67cb25 |
GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef DEBUG
|
|
Packit |
67cb25 |
printf ("got stepc = %g fc = %g\n", stepc, fc);
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Do a line minimisation in the region (xa,fa) (xc,fc) to find an
|
|
Packit |
67cb25 |
intermediate (xb,fb) satisifying fa > fb < fc. Choose an initial
|
|
Packit |
67cb25 |
xb based on parabolic interpolation */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
intermediate_point (fdf, x, p, dir / pnorm, pg,
|
|
Packit |
67cb25 |
stepa, stepc, fa, fc, x1, dx1, gradient, &stepb, &fb;;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (stepb == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_ENOPROG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
minimize (fdf, x, p, dir / pnorm,
|
|
Packit |
67cb25 |
stepa, stepb, stepc, fa, fb, fc, tol,
|
|
Packit |
67cb25 |
x1, dx1, x2, dx, gradient, &(state->step), f, &g1norm);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_memcpy (x, x2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Choose a new conjugate direction for the next step */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->iter = (state->iter + 1) % x->size;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->iter == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_memcpy (p, gradient);
|
|
Packit |
67cb25 |
state->pnorm = g1norm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* p' = g1 - beta * p */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double beta = -pow (g1norm / g0norm, 2.0);
|
|
Packit |
67cb25 |
gsl_blas_dscal (-beta, p);
|
|
Packit |
67cb25 |
gsl_blas_daxpy (1.0, gradient, p);
|
|
Packit |
67cb25 |
state->pnorm = gsl_blas_dnrm2 (p);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->g0norm = g1norm;
|
|
Packit |
67cb25 |
gsl_vector_memcpy (g0, gradient);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef DEBUG
|
|
Packit |
67cb25 |
printf ("updated conjugate directions\n");
|
|
Packit |
67cb25 |
printf ("p: ");
|
|
Packit |
67cb25 |
gsl_vector_fprintf (stdout, p, "%g");
|
|
Packit |
67cb25 |
printf ("g: ");
|
|
Packit |
67cb25 |
gsl_vector_fprintf (stdout, gradient, "%g");
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_multimin_fdfminimizer_type conjugate_fr_type = {
|
|
Packit |
67cb25 |
"conjugate_fr", /* name */
|
|
Packit |
67cb25 |
sizeof (conjugate_fr_state_t),
|
|
Packit |
67cb25 |
&conjugate_fr_alloc,
|
|
Packit |
67cb25 |
&conjugate_fr_set,
|
|
Packit |
67cb25 |
&conjugate_fr_iterate,
|
|
Packit |
67cb25 |
&conjugate_fr_restart,
|
|
Packit |
67cb25 |
&conjugate_fr_free
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_multimin_fdfminimizer_type
|
|
Packit |
67cb25 |
* gsl_multimin_fdfminimizer_conjugate_fr = &conjugate_fr_type;
|