使用OpenMp在内部循环中循环阵列

Question

片段是我尝试完成的通用版本。

program main 
integer, save:: j,Nj,k,Nk,AllocateStatus 
!double precision:: array(100000) 
double precision, dimension(:), allocatable:: array 

INTEGER NTHREADS, TID, OMP_GET_NUM_THREADS, 
+  OMP_GET_THREAD_NUM 

Nj = 100000 
Nk = 100 
allocate(array(Nk),STAT = AllocateStatus) 
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" 

array = 0.0 
!$OMP PARALLEL PRIVATE(NTHREADS, TID) 
!$OMP DO 
do j=1, Nj 

!print *, "id", OMP_GET_THREAD_NUM() 
!DO COMPUTATIONALLY INTENSIVE PART 


    do k=1, Nk 
     array(k)=array(k)+1 
    enddo 

enddo 
!$OMP END DO 
!$OMP END PARALLEL 

print *, array 
stop 
end 

在非的OpenMP版本的array每个元素将是100000使用OpenMP作为片段中，我得到约99000数组元素我不是我需要什么样的变化做出明确获得openmp版本以获得与串行版本相同的输出。

编辑：

在外部循环中完成的操作对相互不依赖的，但这些操作的输出需要在像array一个变量来累积地跟踪。因此，每个线程的临时数组可以在外循环完成后组合起来，但是我不知道如何进行缩减。下面的代码是否有意义？

program main 
integer, save:: j,Nj,k,Nk,AllocateStatus 
!double precision:: array(100000) 
double precision, dimension(:), allocatable:: a_tmp,a_total 

INTEGER NTHREADS, TID, OMP_GET_NUM_THREADS, 
+  OMP_GET_THREAD_NUM 

Nj = 100000 
Nk = 100 
allocate(a_tmp(Nk),STAT = AllocateStatus) 
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" 

    allocate(a_total(Nk),STAT = AllocateStatus) 
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" 

a_tmp = 0.0 
a_total = 0.0 
!$OMP PARALLEL PRIVATE(NTHREADS, TID,a_tmp) 
!$OMP DO 
do j=1, Nj 

!print *, "id", OMP_GET_THREAD_NUM() 

    do k=1, Nk 
     a_tmp(k)=a_tmp(k)+1 
    enddo 

enddo 
!$OMP END DO 
a_total=a_total+a_tmp 
!$OMP END PARALLEL 

print *, a_total 
stop 
end 

Answer 1

这是一个相当普遍的答案。

你的循环很不稳固。正如你写的那样，j的值将在所有线程之间平均分配，所以线程1得到j=1..Nj/num_threads，线程2得到j=(Nj/num_threads)+1..2*Nj/num_threads等等。但是每个线程都会执行

do k=1, Nk 

对于所有值1..Nk。这意味着所有线程将更新array(k)，并且您无法控制何时或如何发生。特别是，你不能阻止线程1读取一个值，然后线程2读取相同的值，然后将它加1并写回，然后线程1为同一元素写回自己的值 - 从而丢失一次更新变量的值。

您已编写了数据竞赛。没有什么可耻的，我们都做到了。

如何删除这是另一个问题，它取决于你想要达到的目标。例如，您可以并行k而不是j。但是您向我们展示的简单代码可能太简单了。我想说，在OpenMP中编写嵌套循环相对比较少见，其中内部循环的内容没有以某种方式包含外部循环迭代器的值（您的案例中为j），因此它们没有实现数据竞赛。

例如

!$OMP DO 
    do j=1, Nj 
     do k=1, Nk 
      array(j)=array(j)+1 
     enddo 
    enddo 

不具有数据争用;编译器/运行时将负责分发循环迭代，并且不会发生多次写入。另一方面，它表达了与您编程不同的操作。

另一种方法是使用OpenMP collapse指令，该指令将融合两个循环并负责消除数据竞争。无论是你的教程资料还是其他Q和As都会告诉你如何使用它。