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:
- Load a dynamic library (DLL) containing the routines required.
- Call the routine from the library.
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
SCALEto be of typeLOGICAL. 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, whereNAG_ROUTINEis 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.
