在 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 函数有以下两个步骤:
    1. 加载函数所需的动态链接库 (DLL)。
    2. 自算法库中调用函数。
    第一个步骤可以采用 R 的函数 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 与其他编译器的使用详情,您可以 参考此处

返回