Re: [Bug-gsl] Possible bug in the spline interpolation routines

[Lemoine, Roberto]

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
address@hidden
Tel: +49 (0) 231/755-5171


-----Ursprüngliche Nachricht-----
Von: Lemoine, Roberto 
Gesendet: Dienstag, 20. August 2013 17:44
An: address@hidden
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