OpenMP 対応の Fortran ルーチンをいくつか呼び出す R のプログラムがあります。2 つの Fortran ルーチンsub_1とsub_2. 1 つ目は R 関数で 2 回呼び出され、2 つ目は 1 回呼び出されます。両方のルーチンは、いくつかの小さな点を除いてほとんど同じです。最初のルーチンを呼び出し、次に 2 番目、次に最初のルーチンを呼び出します。ただし、両方とも openMP を有効にすると、最初の fortran ルーチンを 2 回目に使用すると、関数は何も実行しなくなります (エラーが発生したり、実行が停止したりせず、そのまま待機します)。
openMP を無効にするとsub_1、すべて正常に動作します。代わりに で openMP を無効にするsub_2と、 の 2 回目の使用時に同じ方法で再びハングしますsub_1。これは明らかに最初の使用料を通過するため、奇妙です。
スレッドが適切に閉じていないか何かに関係しているのではないかと思いました（openMPについてあまり知りません）。ただし、別の奇妙な点は、これら 3 つのルーチンを呼び出す R 関数が 4 回呼び出されていることsub_2ですsub_2。なぜこれを行うのか、私にはまったくわかりません！参考までに、これは のコードですsub_1。
subroutine correlation_dd_rad(s_bins,min_s,end_s,n,pos1,dd,r)
!!! INTENT IN !!!!!!!!
integer :: s_bins !Number of separation bins
integer :: N !Number of objects
real(8) :: pos1(3,N) !Cartesian Positions of particles
real(8) :: min_s !The smallest separation calculated.
real(8) :: end_s !The largest separation calculated.
real(8) :: r(N) !The radii of each particle (ascending)
!!! INTENT OUT !!!!!!!
real(8) :: dd(N,s_bins) !The binned data.
!!! LOCAL !!!!!!!!!!!!
integer :: i,j !Iterators
integer :: bin
real(8) :: d !Distance between particles.
real(8) :: dr,mins,ends
real(8),parameter :: pi = 3.14159653589
integer :: counter
dd(:,:) = 0.d0
dr = (end_s-min_s)/s_bins
!Perform the separation binning
mins = min_s**2
ends = end_s**2
counter = 1000
!$OMP parallel do private(d,bin,j)
do i=1,N
!$omp critical (count_it)
counter = counter - 1
!$omp end critical (count_it)
if(counter==0)then
counter = 1000
write(*,*) "Another Thousand"
end if
do j=i+1,N
if(r(j)-r(i) .GT. end_s)then
exit
end if
d=(pos1(1,j)-pos1(1,i))**2+&
&(pos1(2,j)-pos1(2,i))**2+&
&(pos1(3,j)-pos1(3,i))**2
if(d.LT.ends .AND. d.GT.mins)then
d = Sqrt(d)
bin = Floor((d-min_s)/dr)+1
dd(i,bin) = dd(i,bin)+1.d0
dd(j,bin) = dd(j,bin)+1.d0
end if
end do
end do
!$OMP end parallel do
write(*,*) "done"
end subroutine
なぜこれが起こるのか誰にも手がかりがありますか??
乾杯。
問題を再現する、考えられる最小の例を追加します (ちなみに、これは R の問題である必要があります。ここに示したタイプの小さな例ですが、fortran で記述しても問題なく動作します)。したがって、上記のコードと次のコードを fortran で共有オブジェクトにコンパイルしていますcorrelate.so。
subroutine correlation_dr_rad(s_bins,min_s,end_s,n,pos1,n2,pos2,dd,r1,r2)
!!! INTENT IN !!!!!!!!
integer :: s_bins !Number of separation bins
integer :: N !Number of objects
integer :: n2
real(8) :: pos1(3,N) !Cartesian Positions of particles
real(8) :: pos2(3,n2) !random particles
real(8) :: end_s !The largest separation calculated.
real(8) :: min_s !The smallest separation
real(8) :: r1(N),r2(N2) !The radii of particles (ascending)
!!! INTENT OUT !!!!!!!
real(8) :: dd(N,s_bins) !The binned data.
!!! LOCAL !!!!!!!!!!!!
integer :: i,j !Iterators
integer :: bin
real(8) :: d !Distance between particles.
real(8) :: dr,mins,ends
real(8),parameter :: pi = 3.14159653589
integer :: counter
dd(:,:) = 0.d0
dr = (end_s-min_s)/s_bins
!Perform the separation binning
mins = min_s**2
ends = end_s**2
write(*,*) "Got just before parallel dr"
counter = 1000
!$OMP parallel do private(d,bin,j)
do i=1,N
!$OMP critical (count)
counter = counter - 1
!$OMP end critical (count)
if(counter==0)then
write(*,*) "Another thousand"
counter = 1000
end if
do j=1,N2
if(r2(j)-r1(i) .GT. end_s)then
exit
end if
d=(pos1(1,j)-pos2(1,i))**2+&
&(pos1(2,j)-pos2(2,i))**2+&
&(pos1(3,j)-pos2(3,i))**2
if(d.GT.mins .AND. d.LT.ends)then
d = Sqrt(d)
bin = Floor((d-min_s)/dr)+1
dd(i,bin) = dd(i,bin)+1.d0
end if
end do
end do
!$OMP end parallel do
write(*,*) "Done"
end subroutine
次に、R には次の関数があります。最初の 2 つは、上記の fortran コードをラップするだけです。3 番目は、実際のコードと同様の方法で呼び出します。
correlate_dd_rad = function(pos,r,min_r,end_r,bins){
#A wrapper for the fortran routine of the same name.
dyn.load('correlate.so')
out = .Fortran('correlation_dd_rad',
s_bins = as.integer(bins),
min_s = as.double(min_r),
end_s = as.double(end_r),
n = as.integer(length(r)),
pos = as.double(t(pos)),
dd = matrix(0,length(r),bins), #The output matrix.
r = as.double(r))
dyn.unload('correlate.so')
return(out$dd)
}
correlate_dr_rad = function(pos1,r1,pos2,r2,min_r,end_r,bins){
#A wrapper for the fortran routine of the same name
N = length(r1)
N2 = length(r2)
dyn.load('correlate.so')
out = .Fortran('correlation_dr_rad',
s_bins = as.integer(bins),
min_s = as.double(min_r),
end_s = as.double(end_r),
n = N,
pos1 = as.double(t(pos1)),
n2 = N2,
pos2 = as.double(t(pos2)),
dr = matrix(0,nrow=N,ncol=bins),
r1 = as.double(r1),
r2 = as.double(r2))
dyn.unload('correlate.so')
return(out$dr)
}
the_calculation = function(){
#Generate some data to use
pos1 = matrix(rnorm(30000),10000,3)
pos2 = matrix(rnorm(30000),10000,3)
#Find the radii
r1 = sqrt(pos1[,1]^2 + pos1[,2]^2+pos1[,3]^2)
r2 = sqrt(pos2[,1]^2 + pos2[,2]^2+pos2[,3]^2)
#usually sort them but it doesn't matter here.
#Now call the functions
print("Calculating the data-data pairs")
dd = correlate_dd_rad(pos=pos1,r=r1,min_r=0.001,end_r=0.8,bins=15)
print("Calculating the data-random pairs")
dr = correlate_dr_rad(pos1,r1,pos2,r2,min_r=0.001,end_r=0.8,bins=15)
print("Calculating the random-random pairs")
rr = correlate_dd_rad(pos=pos2,r=r2,min_r=0.001,end_r=0.8,bins=15)
#Now we would do something with it but I don't care in this example.
print("Done")
}
これを実行すると、次の出力が得られます。
[1] "Calculating the data-data pairs"
Another Thousand
Another Thousand
Another Thousand
Another Thousand
Another Thousand
Another Thousand
Another Thousand
Another Thousand
Another Thousand
Another Thousand
done
[1] "Calculating the data-random pairs"
Got just before parallel dr
Another thousand
Another thousand
そして、ただそこに座っています... 実際、数回実行すると、ハングする場所が毎回変わることがわかりました。への 2 番目の呼び出しのほとんどを通過する場合もあれば、 へcorrelate_dd_radの呼び出しの途中までしか到達しない場合もありますcorrelate_dr_rad。