OpenMP calling a function gives wrong results

Question

Hi I am trying to put a do loop in different threads. Now inside the do loop I am calling a function which again calls some subroutine and adding to a total sum. Now if I put parallel enclosing the do loop, it is giving random results however I see that if I put the function inside CRITICAL environment it gives the correct result. But this costs more cpu time and does not improve the speed at all. I tested with a small test program and check that my logic is correct. However in a big program (which I can not post here) this only works when I enclose the function call in CRITICAL. Below I give the test program: (my test program works and gives correct result however in the big program I see that funb is not correctly taken in different threads unless it is in CRITICAL environment.)

       sum=0d0
 !$OMP  PARALLEL PRIVATE(i,j,sum1,xcn,fun)
        ithrd=OMP_GET_THREAD_NUM()
 !$OMP DO
       do i=1,5
         sum1=0d0
         do j=1,3
          xcn=i+j+xx
 !$OMP CRITICAL
         fun=funb(xnc)
         write(*,*)fun
 !$OMP END CRITICAL
         sum1=sum1+fun
         enddo
       enddo
 !$OMP END DO
 !$OMP CRITICAL
        sum=sum+sum1
 !$OMP END CRITICAL
 !$OMP END PARALLEL
        write(*,*)sum

If I remove OMP CRITICAL in the big program I see that different threads are taking same values for funb in different threads which should be different. Therefore my understanding is: there is some restriction in the function being called in PARALLEL section. I would be thankful if anybody can clarify the issue.

The function funb given as:

      COMPLEX*16 FUNCTION FUNB(ZAA)

  IMPLICIT COMPLEX*16 (A-H,O-Z)
  real*8 X1,X2
  COMMON/ZVAR/ZA
  COMMON/XVAR/X1,X2
  ZA=ZAA
  call myinvini
  call myinvc(x2,fout)
  funb=fout
  RETURN
  END

myinvini are some data for wl8,xl8 but myinvc is again a subroutine:

  subroutine myinvc(x,f2)
  complex*16 dir,dirc,sta,ss,ssc,cn,cnc,f2,ff,ffc,func
  complex*16 f22,ans
  integer igauss,inte,l,m
  double precision x,range,phi,w,z,zz,zr
  double precision st,st0,zint,xbl,a,b,dli,sli
  double precision cpar,zero
  double precision xl8,wl8,xl32,wl32
  dimension zint(51)
  COMMON/iinte/inte  
  complex*16 cbeta
  common /wgauss/ xl8(8),wl8(8),xl32(32),wl32(32)
  common /ccpar/ cpar

  include 'constants.h'
  igauss = 8
  zero=0.0d0
  range=201.0d0
  phi=3.0d0/4.0d0*pi
  dir=dcmplx(dcos(phi),dsin(phi))
  dirc=dcmplx(dcos(phi),-dsin(phi))
  sta=dcmplx(cpar,zero)
  st =dexp(dlog(range)/dble(inte))
  st0=1.0d0
  zint(1)=zero
  do 11 l=1,inte    
     st0 =st0*st
     zint(l+1)=st0-1.0d0
  11    continue

  ss=dcmplx(zero,zero)
  ssc=dcmplx(zero,zero)
  xbl=dlog(x)

  do 23 l=1,inte  ! inte=5
     a=zint(l)
     b=zint(l+1)
     dli=(b-a)/2.d0
     sli=(b+a)/2.d0

  do 24 m=1,igauss
     if(igauss.eq. 8) w=wl8(m)
     if(igauss.eq.32) w=wl32(m)
     if(igauss.eq. 8) zz=xl8(m)
     if(igauss.eq.32) zz=xl32(m)
     z =dli*zz+sli
     cn=sta+z*dir
     cnc=sta+z*dirc

  ff=func(cn)
  ffc=func(cnc)

  ss=ss+ff*dir*exp(-xbl*cn)*w*dli
  ssc=ssc+ffc*dirc*exp(-xbl*cnc)*w*dli
  24    continue
  23    continue
  f2=(ss+ssc)
  return
  end

Pleaae show us the code of the function. And please about what makes a function thread safe. This was discussed here many times before. — Vladimir F
@VladimirF I have added the function. Is the problem related to definition of private/shared variables? — Boogeyman
You have shared data in common which seems not to be your intent. — tim18

IanH IanH · Accepted Answer · 2017-08-13T23:06:14

In the absence of threadprivate directive, common block variables are shared. The function referenced inside the parallel section modifies such a common block variable, this will cause a data race and is not permitted by the openmp standard.

The code uses implicit typing and implicit specification of the data sharing attributes for most of the variables referenced in the openmp construct. These are apalling from a coding style perspective. The code as shown has one likely variable spelling mistake, which would likely have been avoided if implicit specifications were avoided.

OpenMP calling a function gives wrong results

1 Answers