Annotated Source Code

1/* multimin/conjugate_pr.c
2 * 
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
4 * 
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.
9 * 
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.
14 * 
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.
18 */
19
20/* conjugate_pr.c -- Conjugate gradient Polak-Ribiere algorithm */
21
22/* Modified by Brian Gough to use single iteration structure */
23
24#include <config.h>
25#include <gsl/gsl_multimin.h>
26#include <gsl/gsl_blas.h>
27
28#include "directional_minimize.c"
29
30typedef struct
31{
32  int iter;
33  double step;
34  double max_step;
35  double tol;
36  gsl_vector *x1;
37  gsl_vector *dx1;
38  gsl_vector *x2;
39  double pnorm;
40  gsl_vector *p;
41  double g0norm;
42  gsl_vector *g0;
43}
44conjugate_pr_state_t;
45
46static int
47conjugate_pr_alloc (void *vstate, size_t n)
48{
49  conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
50
51  state->x1 = gsl_vector_calloc (n);
52
53  if (state->x1 == 0)
54    {
55      GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
56    }
57
58  state->dx1 = gsl_vector_calloc (n);
59
60  if (state->dx1 == 0)
61    {
62      gsl_vector_free (state->x1);
63      GSL_ERROR ("failed to allocate space for dx1", GSL_ENOMEM);
64    }
65
66  state->x2 = gsl_vector_calloc (n);
67
68  if (state->x2 == 0)
69    {
70      gsl_vector_free (state->dx1);
71      gsl_vector_free (state->x1);
72      GSL_ERROR ("failed to allocate space for x2", GSL_ENOMEM);
73    }
74
75  state->p = gsl_vector_calloc (n);
76
77  if (state->p == 0)
78    {
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);
83    }
84
85  state->g0 = gsl_vector_calloc (n);
86
87  if (state->g0 == 0)
88    {
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);
94    }
95
96  return GSL_SUCCESS;
97}
98
99static int
100conjugate_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)
103{
104  conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
105
106  state->iter = 0;
107  state->step = step_size;
108  state->max_step = step_size;
109  state->tol = tol;
110
111  GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);
112
113  /* Use the gradient as the initial direction */
114
115  gsl_vector_memcpy (state->p, gradient);
116  gsl_vector_memcpy (state->g0, gradient);
117
118  {
119    double gnorm = gsl_blas_dnrm2 (gradient);
120
121    state->pnorm = gnorm;
122    state->g0norm = gnorm;
123  }
124
125  return GSL_SUCCESS;
126}
127
128static void
129conjugate_pr_free (void *vstate)
130{
131  conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
132
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);
138}
139
140static int
141conjugate_pr_restart (void *vstate)
142{
143  conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
144
145  state->iter = 0;
146  return GSL_SUCCESS;
147}
148
149static int
150conjugate_pr_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
1Start Analysis.
151                      gsl_vector * x, double *f,
152                      gsl_vector * gradient, gsl_vector * dx)
153{
154  conjugate_pr_state_t *state = (conjugate_pr_state_t *) vstate;
155
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;
161
162  double pnorm = state->pnorm;
163  double g0norm = state->g0norm;
164
165  double fa = *f, fb, fc;
166  double dir;
167  double stepa = 0.0, stepb, stepc = state->step, tol = state->tol;
168
169  double g1norm;
170  double pg;
171
172  if (pnorm == 0.0 || g0norm == 0.0)
2Take the false branch.    pnorm == 0.
3Take the false branch.    pnorm == 0. || g0norm == 0.
173    {
174      gsl_vector_set_zero (dx);
175      return GSL_ENOPROG;
176    }
177
178  /* Determine which direction is downhill, +p or -p */
179
180  gsl_blas_ddot (p, gradient, &pg);
4Call a function.    gsl_blas_ddot(p, gradient, &pg)
181
182  dir = (pg >= 0.0) ? +1.0 : -1.0;
12The operand has undefined value.    pg >= 0.
183
184  /* Compute new trial point at x_c= x - step * p, where p is the
185     current direction */
186
187  take_step (x, p, stepc, dir / pnorm, x1, dx);
188
189  /* Evaluate function and gradient at new point xc */
190
191  fc = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);
192
193  if (fc < fa)
194    {
195      /* Success, reduced the function value */
196      state->step = stepc * 2.0;
197      *f = fc;
198      gsl_vector_memcpy (x, x1);
199      GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);
200      return GSL_SUCCESS;
201    }
202
203#ifdef DEBUG
204  printf ("got stepc = %g fc = %g\n", stepc, fc);
205#endif
206
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 */
210
211  intermediate_point (fdf, x, p, dir / pnorm, pg,
212                      stepa, stepc, fa, fc, x1, dx1, gradient, &stepb, &fb);
213
214  if (stepb == 0.0)
215    {
216      return GSL_ENOPROG;
217    }
218
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);
222
223  gsl_vector_memcpy (x, x2);
224
225  /* Choose a new conjugate direction for the next step */
226
227  state->iter = (state->iter + 1) % x->size;
228
229  if (state->iter == 0)
230    {
231      gsl_vector_memcpy (p, gradient);
232      state->pnorm = g1norm;
233    }
234  else
235    {
236      /* p' = g1 - beta * p */
237
238      double g0g1, beta;
239
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) */
243
244      gsl_blas_dscal (-beta, p);
245      gsl_blas_daxpy (1.0, gradient, p);
246      state->pnorm = gsl_blas_dnrm2 (p);
247    }
248
249  state->g0norm = g1norm;
250  gsl_vector_memcpy (g0, gradient);
251
252#ifdef DEBUG
253  printf ("updated conjugate directions\n");
254  printf ("p: ");
255  gsl_vector_fprintf (stdout, p, "%g");
256  printf ("g: ");
257  gsl_vector_fprintf (stdout, gradient, "%g");
258#endif
259
260  return GSL_SUCCESS;
261}
262
263
264
265static const gsl_multimin_fdfminimizer_type conjugate_pr_type = {
266  "conjugate_pr",               /* name */
267  sizeof (conjugate_pr_state_t),
268  &conjugate_pr_alloc,
269  &conjugate_pr_set,
270  &conjugate_pr_iterate,
271  &conjugate_pr_restart,
272  &conjugate_pr_free
273};
274
275const gsl_multimin_fdfminimizer_type
276  * gsl_multimin_fdfminimizer_conjugate_pr = &conjugate_pr_type;
277
Events list
Event 1
Event 2
Event 3
Event 4
Event 12