is romberg's numerical integration thread safe?

General OpenMP discussion

is romberg's numerical integration thread safe?

Postby enrico » Tue Jul 22, 2008 2:22 am

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
enrico
 
Posts: 8
Joined: Wed Jul 16, 2008 3:32 pm

Re: is romberg's numerical integration thread safe?

Postby ejd » Tue Jul 22, 2008 7:13 am

I need a little more information. What is it that you have parallelized? I don't see any OpenMP directives in this code.
ejd
 
Posts: 1025
Joined: Wed Jan 16, 2008 7:21 am

Re: is romberg's numerical integration thread safe?

Postby enrico » Tue Jul 22, 2008 3:23 pm

Hi ejd,

thanks for helping. The code is long, but a strawman may be the following driver which calls the numerical integrator in parallel. I thought it should be straightforward since it is an embarrassing parallelism, but apparently it's not... Sorry, probably the mistake is a silly one...

program test
use integrator
use nrtype
use omp_lib
implicit none
real(dp) :: a, res
integer(I4B) :: n

open(unit=1, file='out.txt')

!$omp parallel do private(n,a)
do n=1,10

a=real(n,dp)
call integrate(a,res)

!$omp critical
write(1,*) res
!$omp end critical

end do
!$omp end parallel do

end program

module integrator
implicit none
contains

function trial_func(x,pars)
use nrtype
real(dp), dimension(2) :: trial_func
real(dp), intent(in) :: x
real(dp), dimension(:), intent(in) :: pars

trial_func(1)=pars(1)/sqrt(x)
trial_func(2)=0.d0

end function

subroutine integrate(param,res)
use nrtype
use integrals
implicit none
real(dp), intent(in) :: param
real(dp), intent(out) :: res


call romb(trial_func,(/param/), 1.d-3,1.d0,res,1.d-6,*10)

return

10 write(*,*) "aborting..."

stop

end subroutine

end module
enrico
 
Posts: 8
Joined: Wed Jul 16, 2008 3:32 pm

Re: is romberg's numerical integration thread safe?

Postby enrico » Tue Jul 22, 2008 3:38 pm

sorry, the earlier driver is trivially wrong. here's a new one resembling to what I do in the real code

program test
use integrator
use nrtype
use omp_lib
implicit none
real(dp) :: a, res, res0
integer(I4B) :: n

res0=0.d0

!$omp parallel do private(n,a,res)
do n=1,10

a=real(n,dp)
call integrate(a,res)

!$omp critical
res0=res0+res
!$omp end critical

end do
!$omp end parallel do

write(*,*) res0

end program
enrico
 
Posts: 8
Joined: Wed Jul 16, 2008 3:32 pm

Re: is romberg's numerical integration thread safe?

Postby ejd » Wed Jul 23, 2008 6:18 am

Unfortunately I don't see your problem. After I put in the Fortran interface blocks and got the program to compile I tried running it on a Sun Sparc system using Solaris and Sun Studio 12 and saw the following:
Code: Select all
serial version
% a.out
106.52171618919004

OpenMP version
% setenv OMP_DYNAMIC FALSE
% setenv OMP_NUM_THREADS 2
% a.out
106.52171618919007

% setenv OMP_NUM_THREADS 4
% a.out
106.52171618919006

I then tried it using the Intel compiler (9.1.031) on an Intel chipset running RedHat Linux and got the following:
Code: Select all
serial version
% a.out
   106.521716189190

OpenMP version
% setenv OMP_DYNAMIC FALSE
% setenv OMP_NUM_THREADS 2
% a.out
   106.521716189190
% setenv OMP_NUM_THREADS 4
% a.out
   106.521716189190

I don't see the variation you described nor do I see any races or other problems.
ejd
 
Posts: 1025
Joined: Wed Jan 16, 2008 7:21 am

Re: is romberg's numerical integration thread safe?

Postby enrico » Wed Jul 23, 2008 9:01 am

yes, I think you're right. I did the same on my machine and I don't get the problem with this simple code. Sorry about it, but what I do in the full code is essentially that, just in a more involved way. I have an ode integrator routine which uses a derivative routine which in turns uses the integration routines that I showed. No saved variables, explicit interfaces, embarrassing parallelism, but the code is executed differently by different threads. So differently that sometimes the integrals are computed alright, while sometimes the integration fails because the accuracy remains too low. If I eliminate the integrals I seem to have no problems, but that may not be the problem as we've seen... Don't know what to think, I'll probably give up as the time spent on this is probably exceeding the benefit :-). Thanks anyway
enrico
 
Posts: 8
Joined: Wed Jul 16, 2008 3:32 pm


Return to Using OpenMP

Who is online

Users browsing this forum: Google [Bot] and 8 guests