[Bug-gsl] GSL_SIMAN Bug

[Zach Jones]

To whom this may concern,

I am using GSL version 1.5 and have tested the bug with version 1.6.

I am using the function gls_siman_solve with a 8-parameter system.  On 
iteration 1246 I get the following message: "*** glibc detected *** 
free(): invalid next size (fast): 0x080fbaf0 *** \\ Aborted."

I have changed different variables and modified the error function, but 
the problem persists at the same iteration.

Attached is the file responislbe for the call to gsl_siman_solve,  the 
code mentions ErrorFun3,  it is the primary function responsible for 
calculating the error.  It does not use any calls to malloc or free.

Thanks for any help you can provide.

Sincerely,
Zachary H. Jones

#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_siman.h>

/* set up parameters for this simulated annealing run */

/* how many points do we try before stepping */
#define N_TRIES 201             

/* how many iterations for each T? */
#define ITERS_FIXED_T 10        

/* max step size in random walk */
#define STEP_SIZE 1            

/* Boltzmann constant */
#define K 1.0                   

/* initial temperature */
#define T_INITIAL 0.001         

/* damping factor for temperature */
#define MU_T 1.005              
#define T_MIN 2.0e-6

gsl_siman_params_t params 
  = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
     K, T_INITIAL, MU_T, T_MIN};

/* now some functions to test in one dimension */
double ErrorFun3(double*);
double E1(void *xp)
{
//int i;
//   for(i=0; i<8; i++) printf("%f ", ((double*)xp)[i]);
//   printf("\n");
   
   return ErrorFun3((double*)xp);
}

double M1(void *xp, void *yp)
{
   printf("M1\n");
   int i;
   double sqsum=0, dz;
   double *x = (double*)xp;
   double *y = (double*)yp;
   
   for(i=0; i<8; i++) {
      dz = x[i] -y[i];
      sqsum += dz*dz;
   }
      
   return sqrt(sqsum);
}

void S1(const gsl_rng * r, void *xp, double step_size)
{
  double *old_x = ((double *) xp);
   
  int i;
   double u = gsl_rng_uniform(r);
   printf("S1\n", u);
  for(i=0; i<8; i++) {
     old_x[i] = u * 2 * step_size - step_size + old_x[i];
  }

}

void P1(void *xp)
{
  printf ("--%12g", ((double *) xp)[1]);
}
int
siman_start(double* x_initial)
{
  const gsl_rng_type * T;
  gsl_rng * r;
  
  int i;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  gsl_siman_solve(r, x_initial, E1, S1, M1, P1,
                  NULL, NULL, NULL, 
                  sizeof(double), params);
  printf("Finished\n");
  return 0;
}