ifort下OpenMP嵌套DO循环并行化结果异常问题求助

问题核心原因 #

你的代码中，全局变量qs01被所有线程共享，当多个线程同时执行qs01=qs01+...的累加操作时，会出现数据竞争：线程A读取qs01的值还没完成写入，线程B已经读取了旧值并进行计算，最终导致累加结果被覆盖、错乱，这就是线程数增加后结果异常的根本原因。单线程时没有竞争，所以结果正常。

修复方案 #

方法1：使用OpenMP的REDUCTION子句（推荐） #

在PARALLEL DO指令中添加REDUCTION(+:qs01)，让OpenMP自动为每个线程创建qs01的私有副本，并行计算完成后再将所有线程的副本值累加到全局的qs01中，避免竞争。

修改后的代码片段：

C$OMP  PARALLEL DO DEFAULT(SHARED), PRIVATE(kp,i1,kt,j,kf,k,k2), REDUCTION(+:qs01)
      
      do kp=1,np
      do i1=1,nleg1
    do kt=1,nt
      do j=1,nleg
    do kf=1,nf
      do k=1,nleg

      ap=(cp(kp+1)-cp(kp))*dx1(i1)*0.5d0+(cp(kp+1)+cp(kp))*0.5d0
      atp=(ct(kt+1)-ct(kt))*dx2(j)*0.5d0+(ct(kt+1)+ct(kt))*0.5d0
      aphp=(cf(kf+1)-cf(kf))*dx3(k)*0.5d0+(cf(kf+1)+cf(kf))*0.5d0
    

C 'LES COMP DES P ( PX, PY, PZ ) et ( Ki + P )'
      apx=ap*dsin(atp)*dcos(aphp)
      apy=ap*dsin(atp)*dsin(aphp)
      apz=ap*dcos(atp)
    
C VECTEUR Ki+P
      akpx= akix + apx
      akpy= akiy + apy
      akpz= akiz + apz

c       print*,akpx,akpy,akpz
C VECTEUR K01=0.5*(ks-k1) et K01+P


      akp01x=ak01x-apx
      akp01y=ak01y-apy
      akp01z=ak01z-apz
         
C -------------------------------------------------
C--                       ELEMENT DE MATRICE
      t1=tcs(alpha,az,-apx,-apy,-apz,akex,akey,akez,1) *fa(alpha,0) 
      qj01=tcs(al0,az,akpx,akpy,akpz,aksx,aksy,aksz,1)
     $*q0s(al0,az01,akp01x,akp01y,akp01z,ak01x,ak01y,ak01z)

      qj0=q0s(al0,az,akpx,akpy,akpz,aksx,aksy,aksz)
     $*tcs(al0,az01,akp01x,akp01y,akp01z,ak01x,ak01y,ak01z,1)

c   print*,t1
    tj1=t1*(qj01-qj0)

    

C The Integration Over P

    qs01=qs01+aw1(i1)*aw2(j)*aw3(k)*ap*ap*(cf(kf+1)-cf(kf))
     $*(cp(kp+1)-cp(kp))*(ct(kt+1)-ct(kt))*dsin(atp)*tj1/8.d0

      enddo
      enddo
      enddo
    enddo
    enddo
    enddo
C$OMP END PARALLEL DO

方法2：手动使用私有变量累加（兼容旧版编译器） #

如果你的编译器不支持REDUCTION子句，可以手动声明线程私有累加变量，并行计算结束后合并结果：

real*8 qs01, qs01_private ! 声明私有累加变量
C$OMP  PARALLEL DEFAULT(SHARED), PRIVATE(kp,i1,kt,j,kf,k,k2,qs01_private)
      qs01_private = 0.0d0 ! 每个线程初始化私有变量
C$OMP DO
      do kp=1,np
      do i1=1,nleg1
    do kt=1,nt
      do j=1,nleg
    do kf=1,nf
      do k=1,nleg

      ! 中间计算代码不变...

    qs01_private=qs01_private+aw1(i1)*aw2(j)*aw3(k)*ap*ap*(cf(kf+1)-cf(kf))
     $*(cp(kp+1)-cp(kp))*(ct(kt+1)-ct(kt))*dsin(atp)*tj1/8.d0

      enddo
      enddo
      enddo
    enddo
    enddo
    enddo
C$OMP END DO
C$OMP CRITICAL
      qs01 = qs01 + qs01_private ! 临界区安全合并结果
C$OMP END CRITICAL
C$OMP END PARALLEL

额外注意事项 #

• 确认akix, akiy, akiz, ak01x, ak01y, ak01z等变量是只读的共享变量，不会被线程修改，否则也会引发竞争。
• 代码中声明的k2变量未被使用，可以删除该私有声明，避免冗余。

内容的提问来源于stack exchange，提问作者ggt bnd