#include #include #include #include #include #include #include #include #include /* * bug_gnewton.c: Trigger division-by-zero when `f` and `fdf` produce different * results for the same `x`. * Author: Curran D. Muhlberger * Date: 2014-04-27 * Build: `$CC bug_gnewton.c -o bug_gnewton -lgsl -lgslcblas -lm` * Run: `./bug_gnewton` or `GSL_IEEE_MODE=trap-common ./bug_gnewton` */ int my_f(const gsl_vector* x, void* params, gsl_vector* f) { gsl_vector_set(f, 0, gsl_vector_get(x, 0) - 1.0 + GSL_DBL_EPSILON); 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_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(); test_fdfsolver(gsl_multiroot_fdfsolver_gnewton); return EXIT_SUCCESS; }