Octopus
minimizer_low.c
Go to the documentation of this file.
1/*
2 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17 02110-1301, USA.
18
19*/
20
21#include <config.h>
22#include <gsl/gsl_errno.h>
23#include <gsl/gsl_math.h>
24#include <gsl/gsl_min.h>
25#include <gsl/gsl_multimin.h>
26#include <math.h>
27#include <stdio.h>
28#include <stdlib.h>
29
30/* A. First, the interface to the gsl function that calculates the minimum
31 of an N-dimensional function, with the knowledge of the function itself
32 and its gradient. */
33
34/* This is a type used to communicate with Fortran; func_d is the type of the
35 interface to a Fortran subroutine that calculates the function and its
36 gradient. */
37typedef void (*func_d)(const int *, const double *, double *, const int *,
38 double *);
39typedef struct {
42
43/* Compute both f and df together. */
44static void my_fdf(const gsl_vector *v, void *params, double *f,
45 gsl_vector *df) {
46 double *x, *gradient, ff2;
47 int i, dim, getgrad;
48 param_fdf_t *p;
49
50 p = (param_fdf_t *)params;
51
52 dim = v->size;
53 x = (double *)malloc(dim * sizeof(double));
54 gradient = (double *)malloc(dim * sizeof(double));
55
56 for (i = 0; i < dim; i++)
57 x[i] = gsl_vector_get(v, i);
58 getgrad = (df == NULL) ? 0 : 1;
60 p->func(&dim, x, &ff2, &getgrad, gradient);
61
62 if (f != NULL)
63 *f = ff2;
64
65 if (df != NULL) {
66 for (i = 0; i < dim; i++)
67 gsl_vector_set(df, i, gradient[i]);
68 }
69
70 free(x);
71 free(gradient);
72}
74static double my_f(const gsl_vector *v, void *params) {
75 double val;
77 my_fdf(v, params, &val, NULL);
78 return val;
81/* The gradient of f, df = (df/dx, df/dy). */
82static void my_df(const gsl_vector *v, void *params, gsl_vector *df) {
83 my_fdf(v, params, NULL, df);
85
86typedef void (*print_f_ptr)(const int *, const int *, const double *,
87 const double *, const double *, const double *);
88
89int FC_FUNC_(oct_minimize,
90 OCT_MINIMIZE)(const int *method, const int *dim, double *point,
91 const double *step, const double *line_tol,
92 const double *tolgrad, const double *toldr,
93 const int *maxiter, func_d f,
94 const print_f_ptr write_info, double *minimum) {
95 int iter = 0;
96 int status;
97 double maxgrad, maxdr;
98 int i;
99 double *oldpoint;
100 double *grad;
102 const gsl_multimin_fdfminimizer_type *T = NULL;
103 gsl_multimin_fdfminimizer *s;
104 gsl_vector *x;
105 gsl_vector *absgrad, *absdr;
106 gsl_multimin_function_fdf my_func;
107
108 param_fdf_t p;
109
110 p.func = f;
111
112 oldpoint = (double *)malloc(*dim * sizeof(double));
113 grad = (double *)malloc(*dim * sizeof(double));
115 my_func.f = &my_f;
116 my_func.df = &my_df;
117 my_func.fdf = &my_fdf;
118 my_func.n = *dim;
119 my_func.params = (void *)&p;
120
121 /* Starting point */
122 x = gsl_vector_alloc(*dim);
123 for (i = 0; i < *dim; i++)
124 gsl_vector_set(x, i, point[i]);
126 /* Allocate space for the gradient */
127 absgrad = gsl_vector_alloc(*dim);
128 absdr = gsl_vector_alloc(*dim);
130 // GSL recommends line_tol = 0.1;
131 switch (*method) {
132 case 1:
133 T = gsl_multimin_fdfminimizer_steepest_descent;
134 break;
135 case 2:
136 T = gsl_multimin_fdfminimizer_conjugate_fr;
137 break;
138 case 3:
139 T = gsl_multimin_fdfminimizer_conjugate_pr;
140 break;
141 case 4:
142 T = gsl_multimin_fdfminimizer_vector_bfgs;
143 break;
144 case 5:
145 T = gsl_multimin_fdfminimizer_vector_bfgs2;
146 break;
148
149 s = gsl_multimin_fdfminimizer_alloc(T, *dim);
151 gsl_multimin_fdfminimizer_set(s, &my_func, x, *step, *line_tol);
152 do {
153 iter++;
154 for (i = 0; i < *dim; i++)
155 oldpoint[i] = point[i];
156
157 /* Iterate */
158 status = gsl_multimin_fdfminimizer_iterate(s);
160 /* Get current minimum, point and gradient */
161 *minimum = gsl_multimin_fdfminimizer_minimum(s);
162 for (i = 0; i < *dim; i++)
163 point[i] = gsl_vector_get(gsl_multimin_fdfminimizer_x(s), i);
164 for (i = 0; i < *dim; i++)
165 grad[i] = gsl_vector_get(gsl_multimin_fdfminimizer_gradient(s), i);
167 /* Compute convergence criteria */
168 for (i = 0; i < *dim; i++)
169 gsl_vector_set(absdr, i, fabs(point[i] - oldpoint[i]));
170 maxdr = gsl_vector_max(absdr);
171 for (i = 0; i < *dim; i++)
172 gsl_vector_set(absgrad, i, fabs(grad[i]));
173 maxgrad = gsl_vector_max(absgrad);
174
175 /* Print information */
176 write_info(&iter, dim, minimum, &maxdr, &maxgrad, point);
178 /* Store infomation for next iteration */
179 for (i = 0; i < *dim; i++)
180 oldpoint[i] = point[i];
182 if (status)
183 break;
184
185 if ((maxgrad <= *tolgrad) || (maxdr <= *toldr))
186 status = GSL_SUCCESS;
187 else
188 status = GSL_CONTINUE;
189 } while (status == GSL_CONTINUE && iter <= *maxiter);
190
191 if (status == GSL_CONTINUE)
192 status = 1025;
194 gsl_multimin_fdfminimizer_free(s);
195 gsl_vector_free(x);
196 gsl_vector_free(absgrad);
197 gsl_vector_free(absdr);
198
199 free(oldpoint);
200 free(grad);
201
202 return status;
203}
204
205/* B. Second, the interface to the gsl function that calculates the minimum
206 of a one-dimensional function, with the knowledge of the function itself,
207 but not its gradient. */
208
209/* This is a type used to communicate with Fortran; func1 is the type of the
210 interface to a Fortran subroutine that calculates the function and its
211 gradient. */
212typedef void (*func1)(const double *, double *);
213typedef struct {
214 func1 func;
215} param_f1_t;
216
217double fn1(double x, void *params) {
218 param_f1_t *p = (param_f1_t *)params;
219 double fx;
220 p->func(&x, &fx);
221 return fx;
222}
223
224void FC_FUNC_(oct_1dminimize, OCT_1DMINIMIZE)(double *a, double *b, double *m,
225 func1 f, int *status) {
226 int iter = 0;
227 int max_iter = 100;
228 const gsl_min_fminimizer_type *T;
229 gsl_min_fminimizer *s;
230 gsl_function F;
231 param_f1_t p;
232
233 p.func = f;
234
235 F.function = &fn1;
236 F.params = (void *)&p;
237
238 T = gsl_min_fminimizer_brent;
239 s = gsl_min_fminimizer_alloc(T);
240
241 *status = gsl_min_fminimizer_set(s, &F, *m, *a, *b);
242
243 gsl_set_error_handler_off();
245 do {
246 iter++;
247 *status = gsl_min_fminimizer_iterate(s);
248
249 *m = gsl_min_fminimizer_x_minimum(s);
250 *a = gsl_min_fminimizer_x_lower(s);
251 *b = gsl_min_fminimizer_x_upper(s);
252
253 *status = gsl_min_test_interval(*a, *b, 0.00001, 0.0);
254
255 /*if (*status == GSL_SUCCESS) printf ("Converged:\n");*/
256 /*printf ("%5d [%.7f, %.7f] %.7f \n", iter, *a, *b,*m);*/
257 } while (*status == GSL_CONTINUE && iter < max_iter);
258 gsl_min_fminimizer_free(s);
259}
260
261/* C. Third, the interface to the gsl function that calculates the minimum
262 of an N-dimensional function, with the knowledge of the function itself,
263 but not its gradient. */
264
265/* This is a type used to communicate with Fortran; funcn is the type of the
266 interface to a Fortran subroutine that calculates the function and its
267 gradient. */
268typedef void (*funcn)(int *, double *, double *);
269typedef struct {
270 funcn func;
271} param_fn_t;
272
273double fn(const gsl_vector *v, void *params) {
274 double val;
275 double *x;
276 int i, dim;
277 param_fn_t *p;
278
279 p = (param_fn_t *)params;
280 dim = v->size;
281 x = (double *)malloc(dim * sizeof(double));
282
283 for (i = 0; i < dim; i++)
284 x[i] = gsl_vector_get(v, i);
285 p->func(&dim, x, &val);
286
287 free(x);
288 return val;
289}
290
291typedef void (*print_f_fn_ptr)(const int *, const int *, const double *,
292 const double *, const double *);
293
294int FC_FUNC_(oct_minimize_direct,
295 OCT_MINIMIZE_DIRECT)(const int *method, const int *dim,
296 double *point, const double *step,
297 const double *toldr, const int *maxiter,
298 funcn f, const print_f_fn_ptr write_info,
299 double *minimum) {
300 int iter = 0, status, i;
301 double size;
302
303 const gsl_multimin_fminimizer_type *T = NULL;
304 gsl_multimin_fminimizer *s = NULL;
305 gsl_vector *x, *ss;
306 gsl_multimin_function my_func;
307
308 param_fn_t p;
309 p.func = f;
310
311 my_func.f = &fn;
312 my_func.n = *dim;
313 my_func.params = (void *)&p;
314
315 /* Set the initial vertex size vector */
316 ss = gsl_vector_alloc(*dim);
317 gsl_vector_set_all(ss, *step);
318
319 /* Starting point */
320 x = gsl_vector_alloc(*dim);
321 for (i = 0; i < *dim; i++)
322 gsl_vector_set(x, i, point[i]);
323
324 switch (*method) {
325 case 6:
326 T = gsl_multimin_fminimizer_nmsimplex;
327 break;
328 }
329
330 s = gsl_multimin_fminimizer_alloc(T, *dim);
331 gsl_multimin_fminimizer_set(s, &my_func, x, ss);
332
333 do {
334 iter++;
335 status = gsl_multimin_fminimizer_iterate(s);
336
337 if (status)
338 break;
339
340 *minimum = gsl_multimin_fminimizer_minimum(s);
341 for (i = 0; i < *dim; i++)
342 point[i] = gsl_vector_get(gsl_multimin_fminimizer_x(s), i);
343
344 size = gsl_multimin_fminimizer_size(s);
345 status = gsl_multimin_test_size(size, *toldr);
347 write_info(&iter, dim, minimum, &size, point);
348
349 } while (status == GSL_CONTINUE && iter < *maxiter);
350
351 if (status == GSL_CONTINUE)
352 status = 1025;
353
354 gsl_vector_free(x);
355 gsl_vector_free(ss);
356 gsl_multimin_fminimizer_free(s);
357 return status;
358}
real(real64) function func(r1, rn, n, a)
Definition: logrid.F90:355
double fn(const gsl_vector *v, void *params)
static void my_fdf(const gsl_vector *v, void *params, double *f, gsl_vector *df)
double fn1(double x, void *params)
double fabs(double __x) __attribute__((__nothrow__
void(* print_f_fn_ptr)(const int *, const int *, const double *, const double *, const double *)
void(* func_d)(const int *, const double *, double *, const int *, double *)
void(* func1)(const double *, double *)
int FC_FUNC_(oct_minimize, OCT_MINIMIZE) const
static double my_f(const gsl_vector *v, void *params)
void(* funcn)(int *, double *, double *)
static void my_df(const gsl_vector *v, void *params, gsl_vector *df)
void(* print_f_ptr)(const int *, const int *, const double *, const double *, const double *, const double *)
integer, parameter, public minimum
integer, public maxiter
Definition: poisson_cg.F90:141
real(real64) function s()
ptrdiff_t i
Definition: operate_inc.c:12
static double f(double w, void *p)