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.
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.
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
below.
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.
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.
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.
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.
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.
