openmp fortran reduction和critical不适用于数组

Question

我目前正在尝试获得一个fortran FE (有限元)代码来与openmp一起工作。我有一个遍历所有元素的循环，我想要并行工作的ie。以下是代码的一个不起作用的简化部分

     !$omp parallel do default(none) shared(nelm,A,res,enod) private(ie,Fe,B,edof)     
    do ie=1,nelm 
       call calcB(B,A(:,ie))
       call calcFe(Fe,B)
       write(*,*) Fe !writes Fe=40d0, this is correct

        call getEdof(edof,enod(:,ie))
        !$OMP CRITICAL      
        res(edof)=res(edof)+fe
        !$OMP END CRITICAL 
        enddo 
!$omp end parallel do

代码的目的是计算力Fe，然后将其添加到edof的res中。使用calcFe计算力，计算出的力是正确的，但循环后得到的res是错误的。

如果我将calcFe替换为简单的Fe=40d0，然后将其添加到res中，则循环后的结果是正确的

     !$omp parallel do default(none) shared(nelm,A,res,enod) private(ie,Fe,B,edof)


    do ie=1,nelm 
       call calcB(B,A(:,ie))

       Fe=40d0

        call getEdof(edof,enod(:,ie))
        !$OMP CRITICAL      
        res(edof)=res(edof)+fe
        !$OMP END CRITICAL 
        enddo 
!$omp end parallel do

是什么导致了这个错误？在这两种情况下，Fe=40d0都被声明为private，但其中只有一种给出了正确的结果。我可以使用reduction而不是!$ CRITICAL，但它给出了相同的错误。在程序中也使用了几个大的和稀疏的矩阵，但在循环期间是被动的/不使用的。我的主管以前在使用openmp和稀疏矩阵时遇到过问题，并怀疑它们使用的是相同的内存。如果错误不明显，最好使用哪个调试器？我对fortran、openmp和编程都是新手。

我正在使用ifort编译，我的操作系统是ubuntu。

编辑:添加了可以运行的简化代码，尽管此代码可以在代码中运行有两个循环，一个是并行的，一个是串行的，它们应该会产生相同的结果，res和res2

program main

use omp_lib

implicit none

integer :: ie, nelm,enod(4,50*50),edof(12),i,j,k
double precision ::B(12,8),fe(12),A(12,12,2500),res(2601*3),res2(2601*3),finish,start

  !creates  enod
  i=1
  do j=1,50
    ie=j
    do k=1,50
       nelm=k
       enod(:,i)=(/ 51*(nelm-1)+1+ie-1, 51*(nelm-1)+1+ie, 51*(nelm)+1+ie-1, 51*(nelm)+1+ie /)
        i=i+1
    end do
  end do

   A=1d0
   res2=0d0
   nelm=2500   
   start=omp_get_wtime()          
   !$omp parallel do default(none) shared(nelm,A,enod) private(ie,fe,edof,B) reduction(+:res2)
      do ie=1,nelm
            call calcB(B,A(:,:,ie))
          call calcFe(fe,B) !the calculated fe is always 2304
          !can write fe=2304 to get correct result with real code
          call getEdof(edof,enod(:,ie))
          res2(edof)=res2(edof)+fe
      end do
   !$omp end parallel do
   finish=omp_get_wtime()
    write(*,*) 'time: ', finish-start

   res=0d0
   nelm=2500             
   start=omp_get_wtime()
      do ie=1,nelm
            call calcB(B,A(:,:,ie))
          call calcFe(fe,B)
          call getEdof(edof,enod(:,ie))
          res(edof)=res(edof)+fe
      end do
        finish=omp_get_wtime()
    write(*,*) 'time: ', finish-start
    write(*,*) 'difference: ',sum(res2-res) 
    write(*,*) sum(res2)

stop

end program main

subroutine calcB(B,A)
double precision ::B(12,8),A(12,12),C(12)
integer          ::gp
  C=1d0
  do gp=1,8
    B(:,gp)=matmul(A,C)
  end do
end subroutine calcB


subroutine calcFe(fe,B)
double precision  ::fe(12),B(12,8),D(12,12)
integer           ::gp
  fe=0d0
  D=2d0

  do gp=1,8
    fe=fe+matmul(D,B(:,gp))
  end do


end subroutine calcFe


subroutine getEdof(edof,enod)
    implicit none
    integer,intent(in) :: enod(4)
    integer,intent(out):: edof(12)
    edof=0
    edof(1:3)  =(/ enod(1)*3-2, enod(1)*3-1, enod(1)*3 /)
    edof(4:6)  =(/ enod(2)*3-2, enod(2)*3-1, enod(2)*3 /)
    edof(7:9)  =(/ enod(3)*3-2, enod(3)*3-1, enod(3)*3 /)
    edof(10:12)=(/ enod(4)*3-2, enod(4)*3-1, enod(4)*3 /)

end subroutine getedof

和make文件

    FF = ifort -O3 -openmp
OBJ1  = main.f90
ls: $(FORT_OBJS)
    $(FF) -o exec $(OBJ1)

不幸的是，这段代码可以工作，所以我无法复制这个错误。res2和res以串行和并行的方式计算。在我的实际程序中，为了得到一个常量fe，我把所有的值都赋给了1d0。计算出的fe是正确的，如果我在calcFe后面添加一个write(*,*) fe，我就会看到这些值是正确的。然后，我将这些值添加到res2，并将它们与串行res进行比较。因此，它们之间的差异很大，因此不存在数值舍入误差。如果我简单地在我的主程序中声明fe=2304，我会得到正确的答案，即使在使用write时fe已经是2304。

在我的实际程序中，所有子例程都在不同的模块中，因此我需要特别注意吗？另外，在其中一个模块中使用了一些全局变量，它们是只读的，但由于它们没有在子例程中声明，它们不会自动成为私有变量？这应该没有问题，因为我将用于计算fe的所有变量都设置为一个常量，全局变量不直接用于计算fe

Answer 1

解决了这个问题，当我将-openmp添加到我的模块的makefile中时，它开始工作。显然，这些模块需要用-openmp编译，而不仅仅是主文件。