R studio crashing when calling a Fortran module

Question

I am trying to call a module written in Fortran in R studio, but the latter unexpectedly crashes.

I have tried with other Fortran examples (factorial and MC estimation of pi) and these worked well.

The only difference I see between my code that don't work and those examples is that I used a subroutine wrapper for my pure function (calling other pure function) whereas these examples were only relying on subroutines and the fact my function has a vector as input.

Here is the procedure I have followed (or code typed in):

• R CMD SHLIB mypath/myfile.f90 #this compiles without problem; code supplied below.
• dyn.load('ptoh') #This is just the given name of the wrapper subroutine
• .Fortran('ptoh', dimen=as.integer(dimen), p=as.matrix(p), m=as.integer(m), h=integer(1) )

where I arbitrarily set in R: dimen<-3; p<-c(4,6,7); m=3;

Do you have any idea of why is it crashing? Is the way of building my wrapper correct? Is it a problem linked to bind(C, name="ptoh_") ?

I already thank you very much for your help.

Gilles

---

I am working on a Mac running HighSierra with latest version of R. The code in 'myfile.90' is the following :

module hilbert
  implicit none

contains    
  pure function rotate_right(x, d, dimen)
    integer, intent(in) :: x, d, dimen    
    integer :: rotate_right, tmp, mask 

    mask = 2 ** dimen - 1     
    rotate_right = shiftr(x, d)    
    tmp = shiftl(x, dimen - d)    
    rotate_right = iand(ior(rotate_right, tmp), mask)
  end

  pure function rotate_left(x, d, dimen)    
    integer, intent(in) :: x, d, dimen    
    integer :: rotate_left, tmp, mask

    mask = 2 ** dimen - 1        
    rotate_left = shiftl(x, d)    
    tmp = shiftr(x, dimen - d)    
    rotate_left = iand(ior(rotate_left, tmp), mask)
  end    

  pure function gc(i)    
    integer, intent(in) :: i    
    integer :: gc

    gc = ieor(i, shiftr(i, 1))
  end

  pure function entry_point(i)    
    integer, intent(in) :: i    
    integer :: entry_point

    if(i == 0) then    
       entry_point = 0    
    else    
       entry_point = gc(2 * ((i - 1) / 2))
    end if
  end

  pure function exit_point(i, dimen)    
    integer, intent(in) :: i, dimen    
    integer :: exit_point, mask

    mask = 2 ** dimen - 1
    exit_point = ieor(entry_point(mask - i), 2 ** (dimen - 1))
  end

  pure function inverse_gc(g, dimen)    
    integer, intent(in) :: g, dimen    
    integer :: inverse_gc, j

    inverse_gc = g    
    j = 1    
    do while(j < dimen)    
       inverse_gc = ieor(inverse_gc, shiftr(g, j))    
       j = j + 1    
    end do
  end

  pure function intercube_g(i) result(g)    
    integer, intent(in) :: i    
    integer :: g

    g = trailz(ieor(gc(i), gc(i + 1)))
  end

  pure function intracube_d(i, dimen) result(d)    
    integer, intent(in) :: i, dimen    
    integer :: d

    if(i == 0) then    
       d = 0    
    else if(modulo(i, 2) == 0) then    
       d = modulo(intercube_g(i - 1), dimen)    
    else    
       d = modulo(intercube_g(i), dimen)    
    end if    
  end

  pure function transform(e, d, b, dimen) result(t)    
    integer, intent(in) :: e, d, b, dimen    
    integer :: t

    t = rotate_right(ieor(b, e), d + 1, dimen)    
  end

  pure function inverse_transform(e, d, b, dimen) result(t)    
    integer, intent(in) :: e, d, b, dimen    
    integer :: t    

    t = transform(rotate_right(e, d + 1, dimen), dimen - d - 2, b, dimen)
  end

  pure function ptoh(dimen, p, m) result(h)      
    integer, intent(in) :: dimen, p(dimen), m    
    integer :: h, e, d, i, j, l, w    

    h = 0    
    e = 0    
    d = 2    
    do i = m - 1, 0, -1    
       l = 0    
       do j = 1, dimen    
          l = l + 2 ** (j - 1) * ibits(p(j), i, 1)    
       end do    
       l = transform(e, d, l, dimen)    
       w = inverse_gc(l, dimen)    
       e = ieor(e, rotate_left(entry_point(w), d + 1, dimen))    
       d = modulo(d + intracube_d(w, dimen) + 1, dimen)    
       h = ior(shiftl(h, dimen), w)    
    end do    
  end

  subroutine ptoh_R_wrapper(dimen, p, m, h) bind(C, name="ptoh_")
    integer :: dimen, p(dimen), m, h, ptoh
    external ptoh

    h = ptoh(dimen, p, m)
  end
end

Answer 1

You exposed your wrapper subroutine with the name "ptoh_" in your module, but you are calling it with 'ptoh' in R Studio. Actually, 'ptoh' is the name of the wrapped (but also exposed) function.

So, the error may be coming from passing the argument h=integer(1), once that the function is not expecting this argument.

If that is the cause, to fix it you change your calling in R Studio to:

.Fortran('ptoh_', dimen=as.integer(dimen), p=as.matrix(p), m=as.integer(m), h=integer(1))

---

Edit:

Another thing that may be causing the error is that you are overriding the name ptoh inside your wrapper subroutine.

subroutine ptoh_R_wrapper(dimen, p, m, h) bind(C, name="ptoh_")
  integer :: dimen, p(dimen), m, h, ptoh
  external ptoh

  h = ptoh(dimen, p, m)
end

The function ptoh is already accesible inside all the module, including inside the wrapper subroutine (the external ptoh statement should be removed too).

Change the declaration to this:

subroutine ptoh_R_wrapper(dimen, p, m, h) bind(C, name="ptoh_")
  integer :: dimen, p(dimen), m, h

  h = ptoh(dimen, p, m)
end

And check if the error persists, please.

As a last note, as you are making the subroutine interoperable with C, consider using kind parameters from the intrinsic module iso_c_binding in your interoperable arguments (and maybe the value attribute if applicable).

I don't know enough of R to be sure that it is necessary, but it doesn't hurt. The final code would look like this:

subroutine ptoh_R_wrapper(dimen, p, m, h) bind(C, name="ptoh_")
  use, intrinsic :: iso_c_binding, only: c_int
  integer(c_int) :: dimen, p(dimen), m, h

  h = ptoh(dimen, p, m)
end