minimizer_low.c

Go to the documentation of this file.

/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2, or (at your option)
 any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 02110-1301, USA.
 
*/
 
#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_multimin.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
 
#include "string_f.h"
 
/* A. First, the interface to the gsl function that calculates the minimum
   of an N-dimensional function, with the knowledge of the function itself
   and its gradient. */
 
/* This is a type used to communicate with Fortran; func_d is the type of the
   interface to a Fortran subroutine that calculates the function and its
   gradient. */
typedef void (*func_d)(const int *, const double *, double *, const int *,
                       double *);
typedef struct {
  func_d func;
} param_fdf_t;
 
/* Compute both f and df together. */
static void my_fdf(const gsl_vector *v, void *params, double *f,
                   gsl_vector *df) {
  double *x, *gradient, ff2;
  int i, dim, getgrad;
  param_fdf_t *p;
 
  p = (param_fdf_t *)params;
 
  dim = v->size;
  x = (double *)malloc(dim * sizeof(double));
  gradient = (double *)malloc(dim * sizeof(double));
 
  for (i = 0; i < dim; i++)
    x[i] = gsl_vector_get(v, i);
  getgrad = (df == NULL) ? 0 : 1;
 
  p->func(&dim, x, &ff2, &getgrad, gradient);
 
  if (f != NULL)
    *f = ff2;
 
  if (df != NULL) {
    for (i = 0; i < dim; i++)
      gsl_vector_set(df, i, gradient[i]);
  }
 
  free(x);
  free(gradient);
}
 
static double my_f(const gsl_vector *v, void *params) {
  double val;
 
  my_fdf(v, params, &val, NULL);
  return val;
}
 
/* The gradient of f, df = (df/dx, df/dy). */
static void my_df(const gsl_vector *v, void *params, gsl_vector *df) {
  my_fdf(v, params, NULL, df);
}
 
typedef void (*print_f_ptr)(const int *, const int *, const double *,
                            const double *, const double *, const double *);
 
int FC_FUNC_(oct_minimize,
             OCT_MINIMIZE)(const int *method, const int *dim, double *point,
                           const double *step, const double *line_tol,
                           const double *tolgrad, const double *toldr,
                           const int *maxiter, func_d f,
                           const print_f_ptr write_info, double *minimum) {
  int iter = 0;
  int status;
  double maxgrad, maxdr;
  int i;
  double *oldpoint;
  double *grad;
 
  const gsl_multimin_fdfminimizer_type *T = NULL;
  gsl_multimin_fdfminimizer *s;
  gsl_vector *x;
  gsl_vector *absgrad, *absdr;
  gsl_multimin_function_fdf my_func;
 
  param_fdf_t p;
 
  p.func = f;
 
  oldpoint = (double *)malloc(*dim * sizeof(double));
  grad = (double *)malloc(*dim * sizeof(double));
 
  my_func.f = &my_f;
  my_func.df = &my_df;
  my_func.fdf = &my_fdf;
  my_func.n = *dim;
  my_func.params = (void *)&p;
 
  /* Starting point */
  x = gsl_vector_alloc(*dim);
  for (i = 0; i < *dim; i++)
    gsl_vector_set(x, i, point[i]);
 
  /* Allocate space for the gradient */
  absgrad = gsl_vector_alloc(*dim);
  absdr = gsl_vector_alloc(*dim);
 
  // GSL recommends line_tol = 0.1;
  switch (*method) {
  case 1:
    T = gsl_multimin_fdfminimizer_steepest_descent;
    break;
  case 2:
    T = gsl_multimin_fdfminimizer_conjugate_fr;
    break;
  case 3:
    T = gsl_multimin_fdfminimizer_conjugate_pr;
    break;
  case 4:
    T = gsl_multimin_fdfminimizer_vector_bfgs;
    break;
  case 5:
    T = gsl_multimin_fdfminimizer_vector_bfgs2;
    break;
  }
 
  s = gsl_multimin_fdfminimizer_alloc(T, *dim);
 
  gsl_multimin_fdfminimizer_set(s, &my_func, x, *step, *line_tol);
  do {
    iter++;
    for (i = 0; i < *dim; i++)
      oldpoint[i] = point[i];
 
    /* Iterate */
    status = gsl_multimin_fdfminimizer_iterate(s);
 
    /* Get current minimum, point and gradient */
    *minimum = gsl_multimin_fdfminimizer_minimum(s);
    for (i = 0; i < *dim; i++)
      point[i] = gsl_vector_get(gsl_multimin_fdfminimizer_x(s), i);
    for (i = 0; i < *dim; i++)
      grad[i] = gsl_vector_get(gsl_multimin_fdfminimizer_gradient(s), i);
 
    /* Compute convergence criteria */
    for (i = 0; i < *dim; i++)
      gsl_vector_set(absdr, i, fabs(point[i] - oldpoint[i]));
    maxdr = gsl_vector_max(absdr);
    for (i = 0; i < *dim; i++)
      gsl_vector_set(absgrad, i, fabs(grad[i]));
    maxgrad = gsl_vector_max(absgrad);
 
    /* Print information */
    write_info(&iter, dim, minimum, &maxdr, &maxgrad, point);
 
    /* Store infomation for next iteration */
    for (i = 0; i < *dim; i++)
      oldpoint[i] = point[i];
 
    if (status)
      break;
 
    if ((maxgrad <= *tolgrad) || (maxdr <= *toldr))
      status = GSL_SUCCESS;
    else
      status = GSL_CONTINUE;
  } while (status == GSL_CONTINUE && iter <= *maxiter);
 
  if (status == GSL_CONTINUE)
    status = 1025;
 
  gsl_multimin_fdfminimizer_free(s);
  gsl_vector_free(x);
  gsl_vector_free(absgrad);
  gsl_vector_free(absdr);
 
  free(oldpoint);
  free(grad);
 
  return status;
}
 
/* B. Second, the interface to the gsl function that calculates the minimum
   of a one-dimensional function, with the knowledge of the function itself,
   but not its gradient. */
 
/* This is a type used to communicate with Fortran; func1 is the type of the
   interface to a Fortran subroutine that calculates the function and its
   gradient. */
typedef void (*func1)(const double *, double *);
typedef struct {
  func1 func;
} param_f1_t;
 
double fn1(double x, void *params) {
  param_f1_t *p = (param_f1_t *)params;
  double fx;
  p->func(&x, &fx);
  return fx;
}
 
void FC_FUNC_(oct_1dminimize, OCT_1DMINIMIZE)(double *a, double *b, double *m,
                                              func1 f, int *status) {
  int iter = 0;
  int max_iter = 100;
  const gsl_min_fminimizer_type *T;
  gsl_min_fminimizer *s;
  gsl_function F;
  param_f1_t p;
 
  p.func = f;
 
  F.function = &fn1;
  F.params = (void *)&p;
 
  T = gsl_min_fminimizer_brent;
  s = gsl_min_fminimizer_alloc(T);
 
  *status = gsl_min_fminimizer_set(s, &F, *m, *a, *b);
 
  gsl_set_error_handler_off();
 
  do {
    iter++;
    *status = gsl_min_fminimizer_iterate(s);
 
    *m = gsl_min_fminimizer_x_minimum(s);
    *a = gsl_min_fminimizer_x_lower(s);
    *b = gsl_min_fminimizer_x_upper(s);
 
    *status = gsl_min_test_interval(*a, *b, 0.00001, 0.0);
 
    /*if (*status == GSL_SUCCESS) printf ("Converged:\n");*/
    /*printf ("%5d [%.7f, %.7f] %.7f \n", iter, *a, *b,*m);*/
  } while (*status == GSL_CONTINUE && iter < max_iter);
  gsl_min_fminimizer_free(s);
}
 
/* C. Third, the interface to the gsl function that calculates the minimum
   of an N-dimensional function, with the knowledge of the function itself,
   but not its gradient. */
 
/* This is a type used to communicate with Fortran; funcn is the type of the
   interface to a Fortran subroutine that calculates the function and its
   gradient. */
typedef void (*funcn)(int *, double *, double *);
typedef struct {
  funcn func;
} param_fn_t;
 
double fn(const gsl_vector *v, void *params) {
  double val;
  double *x;
  int i, dim;
  param_fn_t *p;
 
  p = (param_fn_t *)params;
  dim = v->size;
  x = (double *)malloc(dim * sizeof(double));
 
  for (i = 0; i < dim; i++)
    x[i] = gsl_vector_get(v, i);
  p->func(&dim, x, &val);
 
  free(x);
  return val;
}
 
typedef void (*print_f_fn_ptr)(const int *, const int *, const double *,
                               const double *, const double *);
 
int FC_FUNC_(oct_minimize_direct,
             OCT_MINIMIZE_DIRECT)(const int *method, const int *dim,
                                  double *point, const double *step,
                                  const double *toldr, const int *maxiter,
                                  funcn f, const print_f_fn_ptr write_info,
                                  double *minimum) {
  int iter = 0, status, i;
  double size;
 
  const gsl_multimin_fminimizer_type *T = NULL;
  gsl_multimin_fminimizer *s = NULL;
  gsl_vector *x, *ss;
  gsl_multimin_function my_func;
 
  param_fn_t p;
  p.func = f;
 
  my_func.f = &fn;
  my_func.n = *dim;
  my_func.params = (void *)&p;
 
  /* Set the initial vertex size vector */
  ss = gsl_vector_alloc(*dim);
  gsl_vector_set_all(ss, *step);
 
  /* Starting point */
  x = gsl_vector_alloc(*dim);
  for (i = 0; i < *dim; i++)
    gsl_vector_set(x, i, point[i]);
 
  switch (*method) {
  case 6:
    T = gsl_multimin_fminimizer_nmsimplex;
    break;
  }
 
  s = gsl_multimin_fminimizer_alloc(T, *dim);
  gsl_multimin_fminimizer_set(s, &my_func, x, ss);
 
  do {
    iter++;
    status = gsl_multimin_fminimizer_iterate(s);
 
    if (status)
      break;
 
    *minimum = gsl_multimin_fminimizer_minimum(s);
    for (i = 0; i < *dim; i++)
      point[i] = gsl_vector_get(gsl_multimin_fminimizer_x(s), i);
 
    size = gsl_multimin_fminimizer_size(s);
    status = gsl_multimin_test_size(size, *toldr);
 
    write_info(&iter, dim, minimum, &size, point);
 
  } while (status == GSL_CONTINUE && iter < *maxiter);
 
  if (status == GSL_CONTINUE)
    status = 1025;
 
  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free(s);
  return status;
}