[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>
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;
}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] odeiv2 and OpenMP,
Andre Brand <=