Hello GSL team,
The GSL implementations of the "Hybrid" family of multidimensional
rootfinders from MINPACK contain a bug where they divide by zero when the
initial guess happens to be the root. This affects the "fsolver" types
"hybrid" and "hybrids" as well as the "fdfsolver" types "hybridj" and
"hybridsj". This results in a floating point exception if trapping is
enabled.
I have attached a small program (bug_hybrid.c), based on the manual's
example code, that triggers this bug when GSL_IEEE_MODE=trap-common. The
current output is:
GSL_IEEE_MODE="trap-common"
type = hybrid
Floating point exception (core dumped)
The expected output is:
GSL_IEEE_MODE="trap-common"
type = hybrid
status = success
root = 1.000000e+00
type = hybrids
status = success
root = 1.000000e+00
type = hybridj
status = success
root = 1.000000e+00
type = hybridsj
status = success
root = 1.000000e+00
I'll note that when I originally discovered this bug, a colleague of mine
claimed that the problem is not present in the original FORTRAN
implementation in MINPACK, but rather was introduced in the GSL port.
Thank you for looking into this, and if you have any questions, just let me
know.
- Curran
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_pow_int.h>
#include <gsl/gsl_vector.h>
/*
* bug_hybrid.c: Trigger division-by-zero when initial guess is root.
* Author: Curran D. Muhlberger <[email protected]>
* Date: 2014-04-27
* Build: `$CC bug_hybrid.c -o bug_hybrid -lgsl -lgslcblas -lm`
* Run: `GSL_IEEE_MODE=trap-common ./bug_hybrid`
*/
int my_f(const gsl_vector* x, void* params, gsl_vector* f) {
gsl_vector_set(f, 0, gsl_vector_get(x, 0) - 1.0);
return GSL_SUCCESS;
}
int my_df(const gsl_vector* x, void* params, gsl_matrix* j) {
gsl_matrix_set(j, 0, 0, 1.0);
return GSL_SUCCESS;
}
int my_fdf(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* j) {
gsl_vector_set(f, 0, gsl_vector_get(x, 0) - 1.0);
gsl_matrix_set(j, 0, 0, 1.0);
return GSL_SUCCESS;
}
int test_fsolver(const gsl_multiroot_fsolver_type* type) {
int status;
int iter = 0;
gsl_multiroot_fsolver* solver;
gsl_multiroot_function func = {&my_f, 1, NULL};
gsl_vector* x;
x = gsl_vector_alloc(1);
gsl_vector_set(x, 0, 1.0);
solver = gsl_multiroot_fsolver_alloc(type, 1);
printf("type = %s\n", gsl_multiroot_fsolver_name(solver));
gsl_multiroot_fsolver_set(solver, &func, x);
do {
++iter;
status = gsl_multiroot_fsolver_iterate(solver);
if (status != GSL_SUCCESS) break;
status = gsl_multiroot_test_residual(
gsl_multiroot_fsolver_f(solver), 1.0e-7);
} while ((status == GSL_CONTINUE) && (iter < 1000));
printf(" status = %s\n", gsl_strerror(status));
printf(" root = %e\n", gsl_vector_get(
gsl_multiroot_fsolver_root(solver), 0));
gsl_multiroot_fsolver_free(solver);
gsl_vector_free(x);
return status;
}
int test_fdfsolver(const gsl_multiroot_fdfsolver_type* type) {
int status;
int iter = 0;
gsl_vector* x;
gsl_multiroot_fdfsolver* solver;
gsl_multiroot_function_fdf func = {&my_f, &my_df, &my_fdf, 1, NULL};
x = gsl_vector_alloc(1);
gsl_vector_set(x, 0, 1.0);
solver = gsl_multiroot_fdfsolver_alloc(type, 1);
printf("type = %s\n", gsl_multiroot_fdfsolver_name(solver));
gsl_multiroot_fdfsolver_set(solver, &func, x);
do {
++iter;
status = gsl_multiroot_fdfsolver_iterate(solver);
if (status != GSL_SUCCESS) break;
status = gsl_multiroot_test_residual(
gsl_multiroot_fdfsolver_f(solver), 1.0e-7);
} while ((status == GSL_CONTINUE) && (iter < 1000));
printf(" status = %s\n", gsl_strerror(status));
printf(" root = %e\n", gsl_vector_get(
gsl_multiroot_fdfsolver_root(solver), 0));
gsl_multiroot_fdfsolver_free(solver);
gsl_vector_free(x);
return status;
}
int main() {
/* export GSL_IEEE_MODE=trap-common */
gsl_ieee_env_setup();
/* Floating point exception */
test_fsolver(gsl_multiroot_fsolver_hybrid);
test_fsolver(gsl_multiroot_fsolver_hybrids);
test_fdfsolver(gsl_multiroot_fdfsolver_hybridj);
test_fdfsolver(gsl_multiroot_fdfsolver_hybridsj);
/* ERROR: approximation to Jacobian has collapsed */
/* test_fsolver(gsl_multiroot_fsolver_broyden); */
return EXIT_SUCCESS;
}
