1 /* multimin/conjugate_pr.c
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
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.
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.
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.
20 /* conjugate_pr.c -- Conjugate gradient Polak-Ribiere algorithm */
22 /* Modified by Brian Gough to use single iteration structure */
25 #include <gsl/gsl_multimin.h>
26 #include <gsl/gsl_blas.h>
28 #include "directional_minimize.c"
47 conjugate_pr_alloc (void *vstate, size_t n)
49 conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
51 state->x1 = gsl_vector_calloc (n);
55 GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
58 state->dx1 = gsl_vector_calloc (n);
62 gsl_vector_free (state->x1);
63 GSL_ERROR ("failed to allocate space for dx1", GSL_ENOMEM);
66 state->x2 = gsl_vector_calloc (n);
70 gsl_vector_free (state->dx1);
71 gsl_vector_free (state->x1);
72 GSL_ERROR ("failed to allocate space for x2", GSL_ENOMEM);
75 state->p = gsl_vector_calloc (n);
79 gsl_vector_free (state->x2);
80 gsl_vector_free (state->dx1);
81 gsl_vector_free (state->x1);
82 GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
85 state->g0 = gsl_vector_calloc (n);
89 gsl_vector_free (state->p);
90 gsl_vector_free (state->x2);
91 gsl_vector_free (state->dx1);
92 gsl_vector_free (state->x1);
93 GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
100 conjugate_pr_set (void *vstate, gsl_multimin_function_fdf * fdf,
101 const gsl_vector * x, double *f, gsl_vector * gradient,
102 double step_size, double tol)
104 conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
107 state->step = step_size;
108 state->max_step = step_size;
111 GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);
113 /* Use the gradient as the initial direction */
115 gsl_vector_memcpy (state->p, gradient);
116 gsl_vector_memcpy (state->g0, gradient);
119 double gnorm = gsl_blas_dnrm2 (gradient);
121 state->pnorm = gnorm;
122 state->g0norm = gnorm;
129 conjugate_pr_free (void *vstate)
131 conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
133 gsl_vector_free (state->g0);
134 gsl_vector_free (state->p);
135 gsl_vector_free (state->x2);
136 gsl_vector_free (state->dx1);
137 gsl_vector_free (state->x1);
141 conjugate_pr_restart (void *vstate)
143 conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
150 conjugate_pr_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
151 gsl_vector * x, double *f,
152 gsl_vector * gradient, gsl_vector * dx)
154 conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
156 gsl_vector *x1 = state->x1;
157 gsl_vector *dx1 = state->dx1;
158 gsl_vector *x2 = state->x2;
159 gsl_vector *p = state->p;
160 gsl_vector *g0 = state->g0;
162 double pnorm = state->pnorm;
163 double g0norm = state->g0norm;
165 double fa = *f, fb, fc;
167 double stepa = 0.0, stepb, stepc = state->step, tol = state->tol;
172 if (pnorm == 0.0 || g0norm == 0.0)
174 gsl_vector_set_zero (dx);
178 /* Determine which direction is downhill, +p or -p */
180 gsl_blas_ddot (p, gradient, &pg);
182 dir = (pg >= 0.0) ? +1.0 : -1.0;
184 /* Compute new trial point at x_c= x - step * p, where p is the
187 take_step (x, p, stepc, dir / pnorm, x1, dx);
189 /* Evaluate function and gradient at new point xc */
191 fc = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);
195 /* Success, reduced the function value */
196 state->step = stepc * 2.0;
198 gsl_vector_memcpy (x, x1);
199 GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);
204 printf ("got stepc = %g fc = %g\n", stepc, fc);
207 /* Do a line minimisation in the region (xa,fa) (xc,fc) to find an
208 intermediate (xb,fb) satisifying fa > fb < fc. Choose an initial
209 xb based on parabolic interpolation */
211 intermediate_point (fdf, x, p, dir / pnorm, pg,
212 stepa, stepc, fa, fc, x1, dx1, gradient, &stepb, &fb);
219 minimize (fdf, x, p, dir / pnorm,
220 stepa, stepb, stepc, fa, fb, fc, tol,
221 x1, dx1, x2, dx, gradient, &(state->step), f, &g1norm);
223 gsl_vector_memcpy (x, x2);
225 /* Choose a new conjugate direction for the next step */
227 state->iter = (state->iter + 1) % x->size;
229 if (state->iter == 0)
231 gsl_vector_memcpy (p, gradient);
232 state->pnorm = g1norm;
236 /* p' = g1 - beta * p */
240 gsl_blas_daxpy (-1.0, gradient, g0); /* g0' = g0 - g1 */
241 gsl_blas_ddot(g0, gradient, &g0g1); /* g1g0 = (g0-g1).g1 */
242 beta = g0g1 / (g0norm*g0norm); /* beta = -((g1 - g0).g1)/(g0.g0) */
244 gsl_blas_dscal (-beta, p);
245 gsl_blas_daxpy (1.0, gradient, p);
246 state->pnorm = gsl_blas_dnrm2 (p);
249 state->g0norm = g1norm;
250 gsl_vector_memcpy (g0, gradient);
253 printf ("updated conjugate directions\n");
255 gsl_vector_fprintf (stdout, p, "%g");
257 gsl_vector_fprintf (stdout, gradient, "%g");
265 static const gsl_multimin_fdfminimizer_type conjugate_pr_type = {
266 "conjugate_pr", /* name */
267 sizeof (conjugate_pr_state_t),
270 &conjugate_pr_iterate,
271 &conjugate_pr_restart,
275 const gsl_multimin_fdfminimizer_type
276 * gsl_multimin_fdfminimizer_conjugate_pr = &conjugate_pr_type;