OpenMP 调用函数会给出错误的结果

Question

您好，我正在尝试在不同的线程中放置一个 do 循环。现在在 do 循环中，我正在调用一个函数，该函数再次调用一些子例程并添加到总和中。现在，如果我将并行封闭 do 循环，它会给出随机结果，但是我看到如果我将函数放在 CRITICAL 环境中，它会给出正确的结果。但这会花费更多的 cpu 时间，并且根本不会提高速度。 我用一个小测试程序进行了测试，并检查我的逻辑是否正确。但是在一个大程序（我不能在这里发布）中，这只有在我将函数调用包含在 CRITICAL 中时才有效。 下面我给出了测试程序：（我的测试程序可以工作并给出正确的结果，但是在大程序中，我发现 funb 没有正确地在不同的线程中使用，除非它处于 CRITICAL 环境中。）

       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

如果我在大程序中删除 OMP CRITICAL，我会看到不同线程在不同线程中为 funb 采用相同的值，这些值应该不同。因此我的理解是：在 PARALLEL 部分调用的函数有一些限制。如果有人能澄清这个问题，我将不胜感激。

函数 funb 给出为：

      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 是 wl8,xl8 的一些数据，但 myinvc 又是一个子程序：

  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

Answer 1

在没有 threadprivate 指令的情况下，公共块变量是共享的。并行部分中引用的函数修改了这样一个公共块变量，这将导致数据竞争，并且是 openmp 标准所不允许的。

该代码对 openmp 构造中引用的大多数变量使用隐式类型和隐式规范的数据共享属性。从编码风格的角度来看，这些都令人震惊。所示代码有一个可能的变量拼写错误，如果避免使用隐式规范，则可能会避免这种错误。