[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] Possible bug in the spline interpolation routines
From: |
Lemoine, Roberto |
Subject: |
Re: [Bug-gsl] Possible bug in the spline interpolation routines |
Date: |
Wed, 21 Aug 2013 12:39:37 +0200 |
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