Multidimensional array operation with R Inline and fortran

Question

I want to integrate Fortran code into my R Package. In order to practice, I am trying to use inline with Fortran code. I can create the cfunction, however, When I try to run it, I get segfault

library(inline)
src = "
integer i, j, k
do k = 1, prok
  do j = 1, ncolv
    do i = 1, nrowv
        emis(i, j,k) = veh(i,j) * lkm(i) * ef(j)*pro(k)
     end do
   end do
end do
"
work <- cfunction(sig = signature(
  nrowv="integer",
  ncolv = "integer",
  prok = "integer",
  veh="numeric",
  lkm = "numeric",
  ef = "numeric",
  pro = "numeric",
  emis = "numeric"), 
  implicit = "none", 
  dim = c("", "", "", "(nrowv,ncolv)",
          "(nrowv)","(ncolv)","(prok)","(nrowv,ncolv,prok)"),
  src, 
  language="F95",
  verbose = F)

work(nrowv = 3, ncolv = 2, prok = 3, veh = 1:6, lkm = 1:3, ef = 1:2, pro = 1:4, emis = as.numeric(1))

gives: Segmentation fault (core dumped) I added the print of the work function to see how is the subroutine.

print(work)
An object of class 'CFunc'
function (nrowv, ncolv, prok, veh, lkm, ef, pro, emis) 
.Primitive(".Fortran")(<pointer: 0x7ffac77db100>, nrowv = as.integer(nrowv), 
    ncolv = as.integer(ncolv), prok = as.integer(prok), veh = as.double(veh), 
    lkm = as.double(lkm), ef = as.double(ef), pro = as.double(pro), 
    emis = as.double(emis))
<environment: 0x55ab99679e70>
code:
  1: 
  2:  SUBROUTINE file77f62087457e ( nrowv, ncolv, prok, veh, lkm, ef, pro, emis )
  3: IMPLICIT none
  4: INTEGER nrowv
  5: INTEGER ncolv
  6: INTEGER prok
  7: DOUBLE PRECISION veh(nrowv,ncolv)
  8: DOUBLE PRECISION lkm(nrowv)
  9: DOUBLE PRECISION ef(ncolv)
 10: DOUBLE PRECISION pro(prok)
 11: DOUBLE PRECISION emis(nrowv,ncolv,prok)
 12: 
 13: integer i, j, k
 14: do k = 1, prok
 15:   do j = 1, ncolv
 16:     do i = 1, nrowv
 17: emis(i, j,k) = veh(i,j) * lkm(i) * ef(j)*pro(k)
 18:       end do
 19:     end do
 20: end do
 21: 
 22: RETURN
 23: END
 24:

I can't help you with your R-issue, but I will note that the triple nested do-loops in the Fortran code are inverted (unless you want to stride all over memory with the potential of blowing out the memory cache). Fortran is a column-major language. — evets, Dec 20 '19 at 03:22
The leftmost index should vary the fast in the loops. `emis(1,1,1), emis(2,1,1), emis(3,1,1)` are stored consecutively in memory locations. `emis(1,1,1), emis(1,1,2), emis(1,1,3)` are stored in memory with a stide of `ncolv * nrowv`. — evets, Dec 20 '19 at 03:49

Ralf Stubner · Accepted Answer · 2019-12-20T17:35:23.643

3

You have to be precise with your in and output variables. Both w.r.t. to data type and size:

nrowv, ncolv and prok are integer, but you supply a numeric.
veh, lkm, ef and pro are numeric, but you supply an integer vector.
pro should be of size prok
the output vector should have size nrowv * ncolv * prok, but you only give it a length one vector

Changing all this, I get some output. I have not verified that it is correct:

> work(nrowv = 3L,
+      ncolv = 2L,
+      prok = 3L,
+      veh = as.numeric(1:6),
+      lkm = as.numeric(1:3),
+      ef = as.numeric(1:2),
+      pro = as.numeric(1:3),
+      emis = numeric(3*2*3))
$nrowv
[1] 3

$ncolv
[1] 2

$prok
[1] 3

$veh
[1] 1 2 3 4 5 6

$lkm
[1] 1 2 3

$ef
[1] 1 2

$pro
[1] 1 2 3

$emis
 [1]   1   4   9   8  20  36   2   8  18  16  40  72   3  12  27  24  60 108

BTW, you might be interested in this answer: https://stackoverflow.com/a/52136973/8416610

edited Dec 20 '19 at 17:35

answered Dec 20 '19 at 16:30

Ralf Stubner

26,263
3
40
75

It helped and thanks for noticing the needed input matrix. However, I still don't get the required output. What I'm missing? – Sergio Dec 20 '19 at 17:19
1

@sergio Actually it does not have to be matrix ... See the reworked answer. – Ralf Stubner Dec 20 '19 at 17:36
That is really cool. Many thanks, And I've seen the link. I might try RCPP + Fortran in the future, but I believe I need firstly to have everything OK with Fortran and then Integrate with RCPP – Sergio Dec 20 '19 at 17:40

Multidimensional array operation with R Inline and fortran

1 Answers1