[Bug-gsl] odeiv2 and OpenMP

[Andre Brand]

odeiv2 and OpenMP

Problem: sequential code works, but parallel code _sometimes_ (not always!  
-- I didn't recognize a pattern) results in:

gsl: driver.c:361: ERROR: integration limits and/or step direction not  
consistent

It seems to be a storage problem. (as far as I know, new, malloc are  
threadsafe, but I tried 'critical', anyway)

In the following code, I used exception handling to illustrate when the  
error occurs. The program doesn't return anything apart from the error  
(from time to time). You maybe have to run it a few times (up to 10 and  
more) to get the error (no further compilation necessary).

Since the error doesn't occur always, I find it very hard to trace. So I'm  
thankful for every hint. (Or bugfix :) )
Compiler: gcc (SUSE Linux) 4.5.0 20100604 [gcc-4_5-branch revision 160292]
gsl version 1.15
CPU: Intel i7


#include <cstddef>
#include <exception>
#include <stdexcept>
#include <iostream>
#include <cstdlib>
extern"C"{
#include "gsl/gsl_errno.h"
#include "gsl/gsl_odeiv2.h"
}
#ifdef _OPENMP
#include <omp.h>
#endif

const int dim=20; //array's dimension

int func (double t,const double* y, double* f, void* params) //defining  
trivial function
{
        for(int i=0;i<dim;i++) f[i]=0.0;
        return GSL_SUCCESS;
}

int main (void)
{
        //a few parameters:
        const int Steps=5;
        const double T=100;
        const double Substeps=10;
        const double eps_abs=1;
        const double eps_rel=0;
        const double h=T/(1.0*Steps);
        const double h_start=h/(1.0*Substeps);
        //parallel section begins. Code works, if this line is omitted
        #pragma omp parallel for default(shared)
        for(int k=0;k<100;k++)
        {
		//since all the variables defined in the following should be private,  
the parallel version should be exactly the same as the sequential one
                double* Lsg=new double[dim];
		gsl_odeiv2_system sys={func,NULL,dim,NULL}; //defining system, omitting  
params and jacobian
		gsl_odeiv2_driver*  
dd=gsl_odeiv2_driver_alloc_y_new(&sys,gsl_odeiv2_step_msadams,h_start,eps_abs,eps_rel);
                
                double t=0;
                try
                {
                        for(int i=0;i<Steps;i++)
                        {       
                                const double t_new=t+h;
                                //applying step
                                if(dd->h ==0) printf("Everything is fine up to 
now");
                                int 
status=gsl_odeiv2_driver_apply(dd,&t,t_new,Lsg);
				if(dd->h ==0) printf("... sometimes, this function sets d->h to zero  
at step i=0 "); //results in Error: gsl: driver.c:361: ERROR: integration  
limits and/or step direction not consistent
                                if(status){     
				printf("h= %f , eps_abs= %f, eps_rel= %f \n", dd->h,  
*(static_cast<double*>(dd->c->state)),*(static_cast<double*>(dd->c->state)+1));
                                throw std::runtime_error("usual error");
                                }
                                //end of step
                        }
                }
		catch(const std::runtime_error& exc){printf("Error appears k = %d at  
core %d and step %f \n", k,omp_get_thread_num(),t);} //error appears  
always at first step i=0, but different k
                gsl_odeiv2_driver_free(dd); //program works, if this line is 
omitted
                delete[] Lsg;
        }
        
        return 0;
}