Calling NAG Fortran Routines from R

Introduction

The NAG library contains a large number of numerical and statistical routines. This document gives a brief introduction on one method of accessing these routines from within the statistical package R (www.r-project.org). This tutorial was written using R version 2.7.0 in the windows environment. The approach adopted should, however, be applicable to other platforms.

Details on an alternative approach for accessing the NAG functionality from within R is given here.

General Notes on Calling Fortran Routines from R

  • Calling Fortran routines from within R involves two steps:
    1. Load a dynamic library (DLL) containing the routines required.
    2. Call the routine from the library.
    The first step is performed using the R function dyn.load(), the second using the R function .Fortran()
  • When calling a Fortran function the following type conversions should be borne in mind:
    R storage mode Fortran type
    Integer INTEGER
    Double DOUBLE PRECISION
    Complex DOUBLE COMPLEX
    Character CHARACTER*255
  • R can only call subroutines, not functions. Therefore to use a Fortran function in R a wrapper must be written to convert the function into a subroutine (see example 1 for details).
  • The function .Fortran() automatically converts the routine name supplied as its first argument into lower case, therefore .Fortran("rs13aaf", ...) and .Fortran("RS13AAF", ...) are both looking for the symbol "rs13aaf" in the DLL. No underscores are attached to the routine name.
  • The function .Fortran() returns all of its arguments, other than the first (i.e. the routine name), as a list.
  • Incorrect variable types can cause unpredictable answers, but may not cause any error / warning messages to be displayed.
  • A detailed description of using Fortran (and C) routines from within R is given at http://www.stats.bris.ac.uk/R/doc/manuals/R-exts.pdf

Notes on Calling the NAG routines from R

  • The NAG DLL is compiled using the stdcall calling convention, where as R uses cdecl. Therefore all NAG routines require a wrapper written in Fortran to be able to call them from R. The Fortran wrapper must be compiled into a DLL that uses the cdecl calling convention. These wrapper functions can then be called from R. Details on creating a DLL is compiler specific, an example of how this is done for the Compaq Visual Fortran compiler can be seen at the foot of the page.

    For more information on DLL calling conventions and details on using R with a variety of compilers see http://www.stats.uwo.ca/faculty/murdoch/software/compilingDLLs
  • When loading a DLL that uses the NAG routines (i.e. the DLL containing the wrappers) into R, the following warning may be observed:

    Warning message:
    DLL attempted to change FPU control word from 8001f to 9001f

    This is a common issue when calling DLLs compiled using a visual studio compiler (as per the NAG library) into R and can be safely ignored. Only one of the NAG routines uses a different FPU control word and this routine returns the control word to its original value before exiting.
  • Some of the NAG routines require different types of input parameters other than those tabulated above, for example C02AGF requires the parameter SCALE to be of type LOGICAL. In these cases the wrapper function may need to change the input type (see example 2 for details).
  • The behaviour of the .Fortran() routine when the Fortran routine it is calling takes character arguments is system specific. It is therefore advisable to avoid using character types by altering the input variable in the wrapper routine (see example 3 for one such case).
  • If the NAG routine being called in the wrapper contains a character argument then it may be necessary to specifically control the way the compiler handles such arguments. Most compilers include (at least one) hidden variable in the argument list for a routine which includes a character argument. This hidden argument holds the length of the character string. Some compilers place this next to the character argument, others place it at the end of the argument list. Most compilers will include options to control what behaviour is used. For example, when using the intel compiler, the declaration
    !DEC$ATTRIBUTES MIXED_STR_LEN_ARG :: NAG_ROUTINE
    
    may be required, where NAG_ROUTINE is replaced by the name of the NAG routine being called.
  • A large number of NAG routines have a parameter called IFAIL. On input this parameter indicates how the routine behaves when encountering an error and on output it returns an error code. When calling the NAG routines from within R this parameter should be set to either -1 or 1 as a value of 0 will cause R to terminate.

Example 1: S13AAF

In this example the NAG routine S13AAF is used to obtain the value of an exponential integral.

Fortran Wrapper

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 Code

The data used corresponds to the example given in section 9 of the NAG documentation for S13AAF.

## Call the routine (S13AAF - Example)
.Fortran("rS13AAF", 
         x=as.double(2), 
         z=as.double(0), 
         ifail=as.integer(-1))

Example 2: C02AGF

In this example the NAG routine C02AGF is used to find all the roots of a real polynomial using a variant of Laguerre's Method.

Fortran Wrapper

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 Code

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)

The data used in this example corresponds to that given in section 9.1 of the NAG documentation for C02AGF.

Example 3: G01NAF

In this example the NAG routine G01NAF is used to compute the cumulants and moments of quadratic forms in Normal variates.

Fortran Wrapper

    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 Code

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

The data used in this example corresponds to that in section 9.1 of the NAG documentation for G01NAF.

Using Compaq Visual Fortran and Other Compilers

When creating a DLL to load into R using the Compaq Visual Fortran Compiler, the following DEC attributes need to be set:

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

Attributes C and REFERENCE ensure that the DLL is using the correct calling convention. The ALIAS statement adds an underscore to the subroutine name before exporting it as the R function .Fortan() assumes that an underscore is present. A similar approach should work with most windows visual Fortran compilers.

Details on using other compilers with R can be found here

The standard R installation comes with its own mechanism for compiling C or Fortran code. A brief description of how to use this mechanism is available here.