1 typedef struct
2 {
3 gsl_function_fdf fdf_linear;
4 gsl_multimin_function_fdf *fdf;
5 /* fixed values */
6 const gsl_vector *x;
7 const gsl_vector *g;
8 const gsl_vector *p;
9
10 /* cached values, for x(alpha) = x + alpha * p */
11 double f_alpha;
12 double df_alpha;
13 gsl_vector *x_alpha;
14 gsl_vector *g_alpha;
15
16 /* cache "keys" */
17 double f_cache_key;
18 double df_cache_key;
19 double x_cache_key;
20 double g_cache_key;
21 }
22 wrapper_t;
23
24 static void
moveto(double alpha,wrapper_t * w)25 moveto (double alpha, wrapper_t * w)
26 {
27 if (alpha == w->x_cache_key) /* using previously cached position */
28 {
29 return;
30 }
31
32 /* set x_alpha = x + alpha * p */
33
34 gsl_vector_memcpy (w->x_alpha, w->x);
35 gsl_blas_daxpy (alpha, w->p, w->x_alpha);
36
37 w->x_cache_key = alpha;
38 }
39
40 static double
slope(wrapper_t * w)41 slope (wrapper_t * w) /* compute gradient . direction */
42 {
43 double df;
44 gsl_blas_ddot (w->g_alpha, w->p, &df);
45 return df;
46 }
47
48 static double
wrap_f(double alpha,void * params)49 wrap_f (double alpha, void *params)
50 {
51 wrapper_t *w = (wrapper_t *) params;
52 if (alpha == w->f_cache_key) /* using previously cached f(alpha) */
53 {
54 return w->f_alpha;
55 }
56
57 moveto (alpha, w);
58
59 w->f_alpha = GSL_MULTIMIN_FN_EVAL_F (w->fdf, w->x_alpha);
60 w->f_cache_key = alpha;
61
62 return w->f_alpha;
63 }
64
65 static double
wrap_df(double alpha,void * params)66 wrap_df (double alpha, void *params)
67 {
68 wrapper_t *w = (wrapper_t *) params;
69 if (alpha == w->df_cache_key) /* using previously cached df(alpha) */
70 {
71 return w->df_alpha;
72 }
73
74 moveto (alpha, w);
75
76 if (alpha != w->g_cache_key)
77 {
78 GSL_MULTIMIN_FN_EVAL_DF (w->fdf, w->x_alpha, w->g_alpha);
79 w->g_cache_key = alpha;
80 }
81
82 w->df_alpha = slope (w);
83 w->df_cache_key = alpha;
84
85 return w->df_alpha;
86 }
87
88 static void
wrap_fdf(double alpha,void * params,double * f,double * df)89 wrap_fdf (double alpha, void *params, double *f, double *df)
90 {
91 wrapper_t *w = (wrapper_t *) params;
92
93 /* Check for previously cached values */
94
95 if (alpha == w->f_cache_key && alpha == w->df_cache_key)
96 {
97 *f = w->f_alpha;
98 *df = w->df_alpha;
99 return;
100 }
101
102 if (alpha == w->f_cache_key || alpha == w->df_cache_key)
103 {
104 *f = wrap_f (alpha, params);
105 *df = wrap_df (alpha, params);
106 return;
107 }
108
109 moveto (alpha, w);
110 GSL_MULTIMIN_FN_EVAL_F_DF (w->fdf, w->x_alpha, &w->f_alpha, w->g_alpha);
111 w->f_cache_key = alpha;
112 w->g_cache_key = alpha;
113
114 w->df_alpha = slope (w);
115 w->df_cache_key = alpha;
116
117 *f = w->f_alpha;
118 *df = w->df_alpha;
119 }
120
121 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)122 prepare_wrapper (wrapper_t * w, gsl_multimin_function_fdf * fdf,
123 const gsl_vector * x, double f, const gsl_vector *g,
124 const gsl_vector * p,
125 gsl_vector * x_alpha, gsl_vector *g_alpha)
126 {
127 w->fdf_linear.f = &wrap_f;
128 w->fdf_linear.df = &wrap_df;
129 w->fdf_linear.fdf = &wrap_fdf;
130 w->fdf_linear.params = (void *)w; /* pointer to "self" */
131
132 w->fdf = fdf;
133
134 w->x = x;
135 w->g = g;
136 w->p = p;
137
138 w->x_alpha = x_alpha;
139 w->g_alpha = g_alpha;
140
141 gsl_vector_memcpy(w->x_alpha, w->x);
142 w->x_cache_key = 0.0;
143
144 w->f_alpha = f;
145 w->f_cache_key = 0.0;
146
147 gsl_vector_memcpy(w->g_alpha, w->g);
148 w->g_cache_key = 0.0;
149
150 w->df_alpha = slope(w);
151 w->df_cache_key = 0.0;
152 }
153
154 static void
update_position(wrapper_t * w,double alpha,gsl_vector * x,double * f,gsl_vector * g)155 update_position (wrapper_t * w, double alpha, gsl_vector *x, double *f, gsl_vector *g)
156 {
157 /* ensure that everything is fully cached */
158 { double f_alpha, df_alpha; wrap_fdf (alpha, w, &f_alpha, &df_alpha); } ;
159
160 *f = w->f_alpha;
161 gsl_vector_memcpy(x, w->x_alpha);
162 gsl_vector_memcpy(g, w->g_alpha);
163 }
164
165 static void
change_direction(wrapper_t * w)166 change_direction (wrapper_t * w)
167 {
168 /* Convert the cache values from the end of the current minimisation
169 to those needed for the start of the next minimisation, alpha=0 */
170
171 /* The new x_alpha for alpha=0 is the current position */
172 gsl_vector_memcpy (w->x_alpha, w->x);
173 w->x_cache_key = 0.0;
174
175 /* The function value does not change */
176 w->f_cache_key = 0.0;
177
178 /* The new g_alpha for alpha=0 is the current gradient at the endpoint */
179 gsl_vector_memcpy (w->g_alpha, w->g);
180 w->g_cache_key = 0.0;
181
182 /* Calculate the slope along the new direction vector, p */
183 w->df_alpha = slope (w);
184 w->df_cache_key = 0.0;
185 }
186