Hey guys,

thanks for your answers to my previous thread. I then found that the code was crashing in parallel because I had the fpe0 and traceback options when compiling (I'm using the linux fortran intel compiler). At some point of the code, a number was getting very small and the option fpe0 was correctly setting it to zero in the serial version, but when compiling with openmp and running it, I was getting an integer divide by zero runtime error. Of course with 128 bit reals the "small number" was still resolvable and the error did not appear. The fix was to compile without fpe0 and traceback.

However, I now have another problem. The code compiles and runs, but it is showing racing conditions (the result is slightly different from the serial version and it changes from run to run). Since mine is an embarrassing parallelism (I want to do call the same routine with different parameters, and the various tasks do no communicate until the end, when I sum their results), I guess my routine is not thread safe. I've passed it with the intel thread checker, but it doesn't complain about anything. I don't have any saved variables/common blocks, and my functions have explicit interfaces (functions are contained in modules, so if I call them with the wrong parameters I get a compilation error). I tried to eliminate piece by piece all the code and the problem *seems* to lie in the numerical integration routine below (romberg's method taken from the numerical recipes book), which each thread should execute multiple times (independently of course). If I eliminate the numerical integrations, I *seem* to have the same results at each run. Do you see any obvious reason why this routine should not be thread safe?

Sorry if I paste a bit too much code, but I've been breaking my head on this for two days now...

Another question: are intel fortan intrinsics thread safe? Eg minloc, maxloc, etc. I guess they should, but they're the only things which may go wrong in the code below, I think.

subroutine romb(func,parameters,a,b,qromb,eps,*)

USE nrtype

USE interpolation

IMPLICIT NONE

REAL(DP), INTENT(IN) :: a,b,eps

REAL(dp), dimension(:), intent(in) :: parameters

REAL(DP), intent(out) :: qromb

INTERFACE

FUNCTION func(x,pars)

USE nrtype

REAL(DP), dimension(2) :: func

REAL(DP), INTENT(IN) :: x

REAL(DP), DIMENSION(:), INTENT(IN) :: pars

END FUNCTION func

END INTERFACE

INTEGER(I4B), PARAMETER :: JMAX=25,JMAXP=JMAX+1,K=5,KM=K-1

REAL(DP), DIMENSION(JMAXP) :: h,s

REAL(DP) :: dqromb

INTEGER(I4B) :: j

h(1)=1.d0

do j=1,JMAX

call trapzd(func,parameters,a,b,s(j),j,*10)

if (j >= K) then

call polint(h(j-KM:j),s(j-KM:j),0.d0,qromb,dqromb)

if (abs(dqromb) <= EPS*abs(qromb)) RETURN

end if

s(j+1)=s(j)

h(j+1)=0.25d0*h(j)

end do

write(*,*) "Too many steps in romb:", abs(dqromb), abs(qromb), abs(dqromb)/abs(qromb)

if( abs(dqromb)/abs(qromb)<1.d-3 ) return

10 return 1

END SUBROUTINE

this uses these other routines:

SUBROUTINE polint(xa,ya,x,y,dy)

USE nrtype; USE nrutil, ONLY : assert_eq,nrerror

IMPLICIT NONE

REAL(DP), DIMENSION(:), INTENT(IN) :: xa,ya

REAL(DP), INTENT(IN) :: x

REAL(DP), INTENT(OUT) :: y,dy

INTEGER(I4B) :: m,n,ns

INTEGER(I4B), dimension(1) :: ns_v

REAL(DP), DIMENSION(size(xa)) :: c,d,den,ho

n=assert_eq(size(xa),size(ya),'polint')

c=ya

d=ya

ho=xa-x

ns_v=minloc(abs(x-xa))

ns=ns_v(1)

y=ya(ns)

ns=ns-1

do m=1,n-1

den(1:n-m)=ho(1:n-m)-ho(1+m:n)

if (any(den(1:n-m) == 0.d0)) &

call nrerror('polint: calculation failure')

den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)

d(1:n-m)=ho(1+m:n)*den(1:n-m)

c(1:n-m)=ho(1:n-m)*den(1:n-m)

if (2*ns < n-m) then

dy=c(ns+1)

else

dy=d(ns)

ns=ns-1

end if

y=y+dy

end do

END SUBROUTINE

SUBROUTINE trapzd(func,parameters,a,b,s,n,*)

use nrtype

implicit none

INTEGER(I4B) :: n,it,j

REAL(dp) :: a,b,s,del,my_sum,tnm,x

REAL(dp), dimension(:) :: parameters

REAL(dp), dimension(2) :: f1, f2

INTERFACE

FUNCTION func(x,pars)

USE nrtype

REAL(DP), dimension(2) :: func

REAL(DP), INTENT(IN) :: x

REAL(DP), DIMENSION(:), INTENT(IN) :: pars

END FUNCTION func

END INTERFACE

intent(in) :: a,b,n,parameters

intent(inout) :: s

if (n==1) then

f1=func(a,parameters)

f2=func(b,parameters)

if(f1(2)/=0.d0.or.f2(2)/=0.d0) return 1

s=0.5d0*(b-a)*(f1(1)+f2(1))

else

it=2**(n-2)

tnm=it

del=(b-a)/tnm

x=a+0.5d0*del

my_sum=0.d0

do j=1,it

f1=func(x,parameters)

if(f1(2)/=0.d0) return 1

my_sum=my_sum+f1(1)

x=x+del

end do

s=0.5d0*(s+(b-a)*my_sum/tnm)

end if

END SUBROUTINE