使用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都會告訴你如何使用它。