How to use Fortran95 code in R

Aim: Write a small test subroutine in Fortran using basic functionality and some calls to the International Mathematics and Statistics Library (IMSL) - compile a dll and call the subroutine from R taking some memory allocation problems into account.

Prerequisites: A working installation of Microsoft Visual Studio and Intel Fortran Compiler v.10 w. IMSL and R 2.6
This link may also be useful.

Save the following code as test.f90 in C:\Fortrancode

!---begin code---

subroutine test(M,N,S1,V,NU,S2,D,INVX)

!DEC$ ATTRIBUTES DLLEXPORT :: test
!Routines necessary for BSKS,
!used in calculation of
!modified Bessel Function of second kind
use BSKS_INT
use UMACH_INT


!We use linear_operators for elegant
!basic linear algebra syntax
use linear_operators


implicit none

!input:
!V vector with values used as input in Bessel fct.
!M number of numbers in vector V
!XDUM matrix with dimensions (M,N)
!D dimension of dummy square matrix MX
!NU fractional order used in calculation of
!value of Bessel fct. modulo 1. Must be positive.
integer, intent(in) :: M,N,D
double precision, intent(in) :: NU, V(M)
double precision :: INVX((D*D)), MX(D,D), XDUM(M,N)


!output:
!S1 sum of elements in XDUM
!S2 values of Bessel fct.
!INVX is both an input and an output vector,
!converted to a square matrix of dim. D
!the inverse is calculated and returned as INVX
double precision, intent(out) :: S1,S2(M)


!some dummy variables
integer :: i,j,k
double precision :: XX, BS(N)


!Set all M*N elements in XDUM to 0.5
XDUM=0.5


!Sum all elements in XDUM and return the result to S1
S1=sum(XDUM)


!Calculation of Bessel fct. values
!corresponding to order NU+N and the
!M different arguments found in V.
!The M results are returned in S2.
do i = 1, M
XX=V(i)
call BSKS (NU, XX, N, BS)
S2(i)=BS(N)
end do


!Transform the INVX vector into a D*D array, store result in MX
!Note use of forall - which is a fast alternative to a
!ordinary for loop.
forall (i =1:D, j=1:D)
MX(i,j)=INVX(i+(j-1)*D)
end forall


!Calculate inverse of MX
MX = .i. MX


!and return result in INVX
forall (i =1:D, j=1:D)
INVX(i+(j-1)*D)=MX(i,j)
end forall


end subroutine test

!---end code---

Run Microsoft Visual Studio as administrator
Choose File -> New Project
Select Intel(R) Fortran -> Library
and click on Dynamic-link library
Then give the project a name (test) and click OK
In the right hand side frame in the Visual Studio main window right-click on Source Files
Choose Add -> Existing Item...
Browse and select the test.f90 file and click OK

Project-> Properties
Select Fortran->Optimization
Set Optimization to Maximize Speed
Select Fortran->Data
Set SEQUENCE Types Obey Alignment Rules to Yes (/align:sequence)
Select Fortran->External Procedures
Set Name Case Interpretation to Lower Case (/names:lowercase)
and Append Underscore to External Names to Yes (\assume:underscore)
Go to Fortran-> Command Line
Write /heap-arrays in Additional Options:
Click OK

Setting the /heap-arrays option circumvent problems when calling the subroutine from R.
(Dummy arrays are allocated using the heap and not the stack.)

Since we use routines from the IMSL library (and maybe IntelMKL) we need a few more adjustments (using 32 bit versions):
Tools->Options...->Intel Fortran->General
Add C:\Program Files\VNI\imsl\fnl600\IA32\include\STATIC
to Includes: (click ... add the path to the list and click OK) and

Add C:\Program Files\VNI\imsl\fnl600\IA32\lib
to Libraries: and Executables:
In Project->Properties->Linker->Input
in the "Additional Dependencies" line, add the following IMSL libraries and Intel MKL library according to your requirement (I highlighted the ones I use):
Static link for 32 bit application: imsl.lib imslsuperlu.lib imslhpc_l.lib imsls_err.lib imslmpistub.lib mkl_c.lib libguide.lib Dynamic link for 32 bit application: imslmkl_dll.lib mkl_c_dll.lib libguide40.lib Static link Intel 64 bit application: imsl.lib imslsuperlu.lib imslscalar.lib mkl_em64t.lib libguide.lib imsls_err.lib imslmpistub.lib Dynamic link Intel 64 bit application: imslmkl_dll.lib mkl_dll.lib libguide40.lib
In Linker ->Command
write /force:multiple in Additional Options:

In the main window change the option above the file editor frame from Debug to Release

Choose Build -> Build test (and wait)
The file test.dll should now be available in
C:\Users\hellmund\Documents\Visual Studio 2005\Projects\test\test\Release
(if your Windows username is hellmund)

In Intel Fortran Compiler version 10 it is not possible to change the default selection of a multi-threaded runtime library to a single-threaded. We thus have to provide som dlls in order to use our test.dll with dyn.load in R.
Copy the files libifcoremd.dll and libmmd.dll from C:\Program Files\Intel\Compiler\Fortran\10.0.025\IA32\Lib to
C:\Fortrancode

I have used the following R-code to call the subroutine.
When writing R-code, I prefer the editor Tinn-R

#begin code

#Load the dll file

dyn.load("C:/Fortrancode/test.dll");

#Is the dll file loaded?
is.loaded("test");

#Set some parameters
M=100000;
V=sqrt(1:M);
D=3;
INVX=matrix(c(1,2,3,4,1,6,7,8,1),ncol=3,byrow=F);
N=5;
S1=0.0;
NU=0.5;
S2=rep(0,M);

#Call the routine
F=.Fortran("test",as.integer(M),as.integer(N),as.double(S1),as.double(V),as.double(NU),as.double(S2),as.integer(D),as.double(INVX))

#See the output
F[[3]]
F[[6]][1:10]

F[[8]]
#Check the inversion of the matrix INVX
solve(INVX)

#end code

Comments

AXM said…
Hi Gunner, excellent details about how to use intel fortran compiler with MS Visual Studio and IMSL. I make my own fortran dll following exactly the way you described here.But R reported warning message "DLL attempted to change FPU control word from 8001f to 9001f
" when I try to load the fortran dll into R and I can not use the function included in the dll, although some people said they can use the function even R report the warning. Have ever encountered some problem or do you know how to solve this problem. Thanks a lot!
Gunnar said…
I can only say these instructions worked, when I used them myself. Are you sure you have made changes to project properties in Release mode and not Debug mode - and that all changes to project properties have been applied?
Otherwise have a look at the more recent post on this blog on Intel Fortran and Visual Studio 2008 - it contains a zip file with pictures of properties settings that might be helpful.

Popular Posts