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 | |
30 | typedef 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 | } |
44 | conjugate_pr_state_t; |
45 | |
46 | static int |
47 | conjugate_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 | |
99 | static int |
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) |
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 | |
128 | static void |
129 | conjugate_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 | |
140 | static int |
141 | conjugate_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 | |
149 | static int |
150 | conjugate_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 | |
265 | static 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 | |
275 | const gsl_multimin_fdfminimizer_type |
276 | * gsl_multimin_fdfminimizer_conjugate_pr = &conjugate_pr_type; |
277 | |