title: OpenMP in Fortran in R - program seems to hang for no reason?
body: <p>I have a program in R which calls a couple of Fortran routines, which are openMP-enabled. There are two Fortran routines <code>sub_1</code> and <code>sub_2</code>. The first one is called twice in an R function, while the second is called once. Both routines are almost identical except for a few minor things. I call the first routine, then the second, then the first again. However, if I have both of them openMP-enabled, the function stops doing anything (doesn't have an error or stop execution, just sits there) when it gets to the second time it uses the first fortran routine. </p>
  
  <p>If I disable the openMP in <code>sub_1</code> then everything runs fine. If I instead disable the openMP in <code>sub_2</code>, then it again hangs in the same fashion on the second usage of <code>sub_1</code>. This is odd because it obviously gets through the first usage fine.</p>
  
  <p>I thought it may be to do with the threads not closing properly or something (I don't know too much about openMP). However, another oddity is that the R function that calls these three routines is being called four times, and if I only enable openMP in <code>sub_2</code>, then this works fine (ie. the second, third etc. call to <code>sub_2</code> doesn't hang). I just have no idea why it would do this! For reference, this is the code for <code>sub_1</code>:</p>
  
  <pre><code>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+&amp;
              &amp;(pos1(2,j)-pos1(2,i))**2+&amp;
              &amp;(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
  </code></pre>
  
  <p>Does anyone have any clue why this would happen??</p>
  
  <p>Cheers.</p>
  
  <p>I'll add in the smallest example that I can think of that does reproduce the problem (by the way, this must be an R problem - a small example of the type that I present here but written in fortran works fine). So I have the above code and the following code in fortran, compiled to the shared object <code>correlate.so</code>:</p>
  
  <pre><code>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+&amp;
              &amp;(pos1(2,j)-pos2(2,i))**2+&amp;
              &amp;(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
  </code></pre>
  
  <p>Then in R, I have the following functions - the first two just wrap the above fortran code. The third calls it in a similar way to my actual code:</p>
  
  <pre><code>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")
  }
  </code></pre>
  
  <p>Running this I get the output:</p>
  
  <pre><code> [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
  </code></pre>
  
  <p>And then it just sits there... Actually, running it a few times has shown that it changes where it hangs each time. Sometimes it gets most of the way through the second call to <code>correlate_dd_rad</code> and others it only gets halfway through the call to <code>correlate_dr_rad</code>.</p>