Dear GSL Team, I have identified the problem. I failed to declare the x_spl and y_spl variables as arrays with as many elements as the spline. Now everything works fine. Sorry for the inconvenience.
Best regards, Roberto Lemoine Technische Universität Dortmund Fakultät Bio- und Chemieingenieurwesen Lehrstuhl für Systemdynamik und Prozessführung G2, Raum 3.29 44221 Dortmund [email protected] Tel: +49 (0) 231/755-5171 -----Ursprüngliche Nachricht----- Von: Lemoine, Roberto Gesendet: Dienstag, 20. August 2013 17:44 An: [email protected] Betreff: Possible bug in the spline interpolation routines Dear GSL Team, I am having a problem with the GSL library (version 1.16). I want to perform an interpolation with cubic splines. However, I do not want to use the dynamic allocation-based routines gsl_interp_accel_alloc and gsl_spline_alloc. Instead, I would like to allocate the interpolator, accelerator and spline structures at compile time. I tried to declare the aforementioned structures (to allocate them at compile time) and give them the same processing as in the aforementioned functions. See the programs below (details about the hardware and the compilation are at the end of this message): ****************************************************************************** ****************************************************************************** #include "../config.h" #include "../string.h" #include <iostream> #include <stdlib.h> #include <stdio.h> #include <math.h> using namespace std; #include "gsl/gsl_errno.h" // GSL: Error handling #include "gsl/gsl_interp.h" // GSL: Interpolators #include "gsl/gsl_spline.h" // GSL: Splines // Accelerator of interpolation searches void gsl_interp_accel_alloc_ct(gsl_interp_accel *a_ptr) { // Initialization of the acceleration object a_ptr->cache = 0; a_ptr->hit_count = 0; a_ptr->miss_count = 0; return; } // ---------------------------------------------------------------- // Interpolation object void gsl_interp_alloc_ct(gsl_interp *interp_ptr, const gsl_interp_type *T, size_t size) { // Initialization of the interpolation object if (T->min_size > size){ cerr << "Insufficient number of points for interpolation " << "type" << endl; } interp_ptr->type = T; interp_ptr->size = size; interp_ptr->state = interp_ptr->type->alloc(size); return; } // ---------------------------------------------------------------- // Cubic splines void gsl_spline_alloc_ct(gsl_spline *spline_ptr, gsl_interp *interp_ptr, double *x, double *y, const gsl_interp_type *T, size_t size) { // Initialization of the interpolation object gsl_interp_alloc_ct(interp_ptr, T, size); // Initialization of the spline object spline_ptr->interp = interp_ptr; spline_ptr->x = x; spline_ptr->y = y; spline_ptr->size = size; return; } // ---------------------------------------------------------------- int main (void) { int i, status; double delta; double x_new[92], y_new[92], x_base[10], y_base[10]; printf ("#m=0,S=2\n"); for (i = 0; i < 10; i++) { x_base[i] = i + 0.5 * sin (i); y_base[i] = i + cos (i * i); printf ("%g %g\n", x_base[i], y_base[i]); } printf ("#m=1,S=0\n"); // HERE ARE MY OWN DECLARATIONS gsl_interp_accel acc; gsl_interp_accel *acc_ptr = &acc; gsl_interp spline_interp; gsl_interp *spline_interp_ptr = &spline_interp; gsl_spline spline; gsl_spline *spline_ptr = &spline; double x_spl, y_spl; // HERE ARE MY ROUTINES FOR PROCESSING THE STRUCTURES gsl_interp_accel_alloc_ct(acc_ptr); gsl_spline_alloc_ct(spline_ptr, spline_interp_ptr, &x_spl, &y_spl, gsl_interp_cspline, 10); gsl_spline_init (spline_ptr, x_base, y_base, 10); delta = 0; for (i = 0; i < 92; i++) { x_new[i] = x_base[0] + delta; y_new[i] = gsl_spline_eval (spline_ptr, x_new[i], acc_ptr); printf ("%g %g\n", x_new[i], y_new[i]); delta += 0.1; } return 0; } ****************************************************************************** ****************************************************************************** So far so good. The problem is when the function gsl_spline_init is executed. After exiting the function, the pointer to my spline spline_ptr is not anymore equal to the address of the spline (&spline). I have tracked the function with gdb, and the change in the value of this pointer occurs in the function gsl_spline_init, after calling: memcpy (spline-y, y_array, size * sizeof(double)); After leaving gsl_spline_init, the pointer spline_ptr is not anymore the address to the spline (&spline), and if I continue the execution of the code, I get a segmentation failure. This does not happen if I use the traditional gsl_interp_accel_alloc and gsl_spline_alloc routines (see below an alternative version of the main function where this failure does not happen). ****************************************************************************** ****************************************************************************** int main (void) { int i, status; double delta; double x_new[92], y_new[92], x_base[10], y_base[10]; printf ("#m=0,S=2\n"); for (i = 0; i < 10; i++) { x_base[i] = i + 0.5 * sin (i); y_base[i] = i + cos (i * i); printf ("%g %g\n", x_base[i], y_base[i]); } printf ("#m=1,S=0\n"); gsl_interp_accel *acc_ptr = gsl_interp_accel_alloc (); gsl_spline *spline_ptr = gsl_spline_alloc (gsl_interp_cspline, 10); gsl_spline_init (spline_ptr, x_base, y_base, 10); delta = 0; for (i = 0; i < 92; i++) { x_new[i] = x_base[0] + delta; y_new[i] = gsl_spline_eval (spline_ptr, x_new[i], acc_ptr); printf ("%g %g\n", x_new[i], y_new[i]); delta += 0.1; } gsl_spline_free (spline_ptr); gsl_interp_accel_free (acc_ptr); return 0; } ****************************************************************************** ****************************************************************************** What could be my failure? Or could this be a bug? I really appreciate your help with this issue. Thanks a lot!! Roberto Lemoine Technical University of Dortmund ****************************************************************************** Additional details: - Compiler: gcc 4.4 - Compilation: gcc -g -I${HOME}/user_install/gsl-1.16/include -c spline_test_2_nodynalloc.cpp gcc -L${HOME}/user_install/gsl-1.16/lib spline_test_2_nodynalloc.o -lgsl -lgslcblas -lm -lstdc++ -o spline_test_2_nodynalloc - System: Ubuntu 12.04 in a VMWare Virtual Machine (Host: Windows 7) - Hardware: Dell Optiplex 980 with, 16 GB RAM, Processor Intel i5
