在 R 中调用 NAG Fortran 函数
索引
简介自 R 中调用 Fortran 函数注意事项
自 R 中调用 NAG 函数注意事项
示例一：S13AAF
示例二：C02AGF
示例三：G01NAF
使用 Compaq Visual Fortran 及其他编译器
回呼 (Callback) 函数
简介
NAG 算法库包含了许多的数值与统计函数。本文将介绍如何在 R (www.r-project.org) 的环境中调用这些函数的方法。 文章中是以 Windows 环境中 R 2.7.0 版本作为说明。我们相信此本文所提的方法将也能够应用在不同的版本与环境中。另一个调用 NAG 函数方法的详细步骤请参考 完整说明。
自 R 中调用 Fortran 函数注意事项
- 在 R 中调用 Fortran 函数有以下两个步骤：
- 加载函数所需的动态链接库 (DLL)。
- 自算法库中调用函数。
dyn.load()，第二个步骤则是使用 R 的函数.Fortran() - 当您要调用 Fortran 函数时，您必须要特别注意数据型别的转换：
R 单元格式 Fortran 型别 Integer INTEGER Double DOUBLE PRECISION Complex DOUBLE COMPLEX Character CHARACTER*255 - R 仅能调用子程序 (subroutines)，而非函数 (functions)。因此，在 R 中若要调用 Fortran 的函数，则必须先要将函数包裹在子程序中。(示例一有详细的说明)。
- 函数 .Fortran() 会自动的将第一个参数名称转换成小写，因此
.Fortran("rs13aaf", ...)与.Fortran("RS13AAF", ...)两者都会参考到 DLL 中的 "rs13aaf"。不需要对函数名称加上底线。 - 函数
.Fortran()会回传除了第一个参数 (函数名称) 外的所有参数。 - 不正确的变量型别可能导致无法预期的答案，且不一定会引起 error / warning 的错误讯息。
- 若您要了解更详细的说明，可以参考 http://www.stats.bris.ac.uk/R/doc/manuals/R-exts.pdf
自 R 中调用 NAG 函数注意事项
- NAG 的动态链接库 (DLL) 是以 stdcall 调用协议所编译的，而 R 采用的是 cdecl。因此所有的 NAG 函数都必须透过 Fortran 写一个包裹程序才能在 R 的环境中使用。
此 Fortran 包裹程序必须以 cdecl 调用协议编译成 DLL。如此便可以使用 R 来调用了。如何编译成 DLL 是编译器的工作，您可以参考 Compaq Visual Fortran 编译器的
例子。
关于更多 R 所支持的编译器及其调用协议等讨论，您可以参考 http://www.stats.uwo.ca/faculty/murdoch/software/compilingDLLs - 当在 R 中载入 DLL 时 (指的是包裹后的 Fortran DLL)，您或许会看到以下的警告讯息：
Warning message:
DLL attempted to change FPU control word from 8001f to 9001f
这是在 R 中调用由 Visual Studio 编译器所编译的 DLL 时普遍会显示的讯息，您可以忽略它。 - 有些 NAG 函数需要有别于以上表格所列的输入参数数据型别，例如 C02AGF 需要的参数
SCALE就是LOGICAL型别。 在这些案例中，包裹程序就必须改变输入的型别 (参考示例二)。 - 当使用
.Fortran()函数调用带有字符参数的函数时，其会被视为系统命令。因此我们需要避免在包裹程序中使用字符参数 (参考示例三)。 - 如果在包裹程序中调用的 NAG 函数带有字符参数，则需要特别以编译器来进行控制。
大部分的编译器在函数的参数列中都包含隐藏 (至少一个) 的变量。这个隐藏的参数记载了字符串的长度。有些编译器会将其放在字符参数后，有些则会放在参数列后。大部分的编译器提供选项来控制这样的机制。例如，当使用 Intel 编译器时，声明
!DEC$ATTRIBUTES MIXED_STR_LEN_ARG :: NAG_ROUTINE
其中NAG_ROUTINE是将被调用使用的 NAG 函数。 - 有很多的 NAG 函数有
IFAIL参数。作为输入时，此参数会可用来控制当遇到错误发生时该如何处理。若在 R 中调用 NAG 函数时，这个参数必须被设为 -1 或 1， 若设为 0 则会导致 R 发生错误。
示例一：S13AAF
本示例中，调用 NAG S13AAF 函数求指数积分。Fortran 包裹程序
C Sample wrapper for a function (S13AAF)
SUBROUTINE rS13AAF(X, Z, IFAIL)
C X and IFAIL are the parameters specified in the documentation
C for NAG routine S13AAF. Z is a new output parameter that will
C hold the results of integral.
INTEGER X, IFAIL
DOUBLE PRECISION Z
DOUBLE PRECISION S13AAF
EXTERNAL S13AAF
Z = S13AAF(X, IFAIL)
END
R 程序
输入数据与 NAG 手册中第 9 节的示例相同 S13AAF。
## Call the routine (S13AAF - Example)
.Fortran("rS13AAF",
x=as.double(2),
z=as.double(0),
ifail=as.integer(-1))
示例二：C02AGF
本示例中，调用 NAG C02AGF 函数的求多项式根。Fortran 包裹程序
C Sample wrapper for a routine that expects logical input C02AGF)
SUBROUTINE rC02AGF(A, N, ISCALE, Z, W, IFAIL)
C A, N, Z, W and IFAIL are the parameters specified in the
C documentation for NAG routine C02AGF. ISCALE is an integer
C representation of SCALE
INTEGER ISCALE, N, IFAIL
DOUBLE PRECISION A(N+1), Z(2,N), W(2*(N+1))
LOGICAL SCALE
EXTERNAL C02AGF
IF (ISCALE.EQ.0) THEN
SCALE = .FALSE.
ELSE
SCALE = .TRUE.
END IF
CALL C02AGF(A, N, SCALE, Z, W, IFAIL)
END
R 程序
rC02AGF = function(a, scale=T, format=T) {
## Wrapper for Fortran function rC02AGF
## a - Corresponds to NAG parameter A
## scale - Corresponds to NAG parameter SCALE
## format - Logical indicating whether the roots are formatted for
## output.
n = length(a) - 1
ans = .Fortran("rC02AGF",
a=as.double(a),
n=as.integer(n),
iscale=as.integer(scale),
z=matrix(as.double(0), 2, n),
w=numeric(20*(n+1)),
ifail=as.integer(-1))
if (format) {
## Format the results, if required
ff = apply(ans$z, 1, function(a) {
ans = c(paste("z = ", a[1]), T)
if (a[2] != 0)
ans = c(paste(ans[1], " + / - ", abs(a[2]), "*i"),
(a[2] < 0))
return(ans)
})
return(ff[1,][as.logical(ff[2,])])
}
else
return(ans$z)
}
## Calling the function
rC02AGF(1:6)
输入数据与 NAG 手册中 9.1 节的示例相同 C02AGF。
示例三：G01NAF
本示例中，调用 NAG G01NAF 函数计算统计值。Fortran 包裹程序
SUBROUTINE rG01NAF(IMOM, IMEAN, N, A, LDA, EMU, SIGMA, LDSIG,
+ L, RKUM, RMOM, WK, IFAIL)
C All variables, with the exception of IMOM and IMEAN are the
C same as specified in the documentation for NAG routine G01NAF.
C IMOM and IMEAN are numeric representations of MOM and MEAN
INTEGER N, LDA, LDSIG, L, IFAIL, IMOM, IMEAN
DOUBLE PRECISION A(LDA,N), EMU(*), SIGMA(LDSIG,N)
DOUBLE PRECISION RKUM(L), RMOM(*), WK(3*N*(N+1)/2+N)
CHARACTER*1 MOM, MEAN
EXTERNAL G01NAF
IF (IMOM.EQ.1) THEN
MOM = 'C'
ELSE
MOM = 'M'
END IF
IF (IMEAN.EQ.1) THEN
MEAN = 'Z'
ELSE
MEAN = 'M'
ENDIF
CALL G01NAF(MOM, MEAN, N, A, LDA, EMU, SIGMA, LDSIG, L,
+ RKUM, RMOM, WK, IFAIL)
END
R 程序
rG01NAF = function(mom='C', mean='Z', a, emu, sigma, l, ifail=-1) {
## Wrapper for fortran routine rG01NAF
## mom - Corresponds to NAG parameter MOM
## mean - Corresponds to MAG parameter MEAN
## a - Corresponds to NAG parameter A
## emu - Corresponds to NAG parameter EMU
## sigma - Corresponds to NAG parameter SIGMA
## l - Corresponds to NAG parameter L
## ifail - Corresponds to NAG parameter IFAIL
n = ncol(a)
## Call the fortran routine
ans = .Fortran("rG01NAF",
imom=as.integer(ifelse(mom=='C', 1, 2)),
imean=as.integer(ifelse(mean=='Z', 1, 2)),
n=as.integer(n),
a=matrix(as.double(a), ncol=ncol(a)),
lda=as.integer(nrow(a)),
emu=as.double(emu),
sigma=matrix(as.double(sigma), ncol=ncol(sigma)),
ldsig=as.integer(nrow(sigma)),
l=as.integer(l),
rkum=numeric(l),
rmom=numeric(l),
wk=numeric(3*n*(n+1)/2+n),
ifail=as.integer(ifail))
## Only really interested in the cumulants and moments
aa = cbind(ans$rkum, ans$rmom)
colnames(aa) = c("Cumulants", "Moments")
return(aa)
}
## Generate the data
beta = 0.8
con = 1
n = 10
l = 4
a = matrix(0, n, n)
a[matrix(c(2:n, 1:(n-1)), ncol=2)] = 0.5
emu = numeric(n)
emu[1] = con*beta
sigma = matrix(0, n, n)
sigma[1,1] = 1
for (i in 1:n) {
if (i != 1) {
emu[i] = beta * emu[i-1]
sigma[i, i] = beta * beta * sigma[i-1, i-1] + 1
}
if (i != n) {
for (j in (i+1):n)
sigma[j, i] = beta*sigma[j-1, i]
}
}
## Call the function and output the results
print(paste("N =", ans$n, "BETA =", beta, "CON =", con))
print(rG01NAF('M', 'M', a, emu, sigma, l))
输入数据与 NAG 手册中 9.1 节的示例相同 G01NAF。
使用 Compaq Visual Fortran 及其他编译器
当使用 Compaq Visual Fortran 编译器编译 DLL 给 R 使用时，您必须要设定以下的 DEC 属性：
SUBROUTINE TEST(A)
C The two DEC lines are compiler specific and used to create a C
DLL with the correct calling convention.
!DEC$ ATTRIBUTES DLLEXPORT :: TEST
!DEC$ ATTRIBUTES C, REFERENCE, ALIAS:'TEST_' :: TEST
INTEGER A
END
C 与 REFERENCE 属性能够确保 DLL 将会使用正确的调用协议。ALIAS 叙述在子程序名称后加上了下底线，能够确保在使用
R .Fortan() 函数时都确实呈现。类似的方法应当都能在大部分的 Windows 上的 Fortran 编译器中使用。
关于 R 与其他编译器的使用详情，您可以 参考此处。