使用 R 内联和 fortran 的多维数组操作

Question

我想将Fortran代码集成到我的R包中。为了练习，我正在尝试将inline与Fortran代码一起使用。我可以创建 cfunction，但是，当我尝试运行它时，我得到段错误

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))

给： 分段故障(核心转储( 我添加了工作函数的打印以查看子例程如何。

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: 

Answer 1

你必须精确地使用你的输入和输出变量。数据类型和大小都包含：

• nrowv、ncolv和prok是整数，但您提供的数字。
• veh、lkm、ef和pro是数值，但你提供一个整数向量。
• pro的大小应为prok
• 输出向量的大小应该nrowv * ncolv * prok，但你只给它一个长度一个向量

改变这一切，我得到了一些输出。我尚未验证它是否正确：

> 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

顺便说一句，您可能对这个答案感兴趣：https://stackoverflow.com/a/52136973/8416610