[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] odeiv2 and OpenMP

From: Andre Brand
Subject: [Bug-gsl] odeiv2 and OpenMP
Date: Wed, 05 Sep 2012 15:00:30 +0200
User-agent: Opera Mail/10.63 (Linux)

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>
#include "gsl/gsl_errno.h"
#include "gsl/gsl_odeiv2.h"
#ifdef _OPENMP
#include <omp.h>

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;
                        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 
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
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 
                delete[] Lsg;
        return 0;

reply via email to

[Prev in Thread] Current Thread [Next in Thread]