### 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
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
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
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
to Includes: (click ... add the path to the list and click OK) and

to Libraries: and Executables:
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 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

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

#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

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!
Unknown 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.

### Comorbidity indexes in SQL

Generating Elixhauser comorbidity index from Danish National Health Register as relational database. ( ICD 10 Coding  in SAS) A lookup-table based version of Charlson comorbidity index I made in SQL. A similar approach can be applied to Elixhauser. SELECT V_CPR, MAX(EI1)+MAX(EI2)+MAX(EI3)+MAX(EI4)+MAX(EI5)+ MAX(EI6)+MAX(EI7)+MAX(EI8)+MAX(EI9)+MAX(EI10)+ MAX(EI11)+MAX(EI12)+MAX(EI13)+MAX(EI14)+MAX(EI15)+ MAX(EI16)+MAX(EI17)+MAX(EI18)+MAX(EI19)+MAX(EI20)+ MAX(EI21)+MAX(EI22)+MAX(EI23)+MAX(EI24)+MAX(EI25)+ MAX(EI26)+MAX(EI27)+MAX(EI28)+MAX(EI29)+MAX(EI30)+MAX(EI31) AS Elixhauser FROM (SELECT V_CPR, -- Congestive Heart Failure CASE WHEN DIAG LIKE 'DI099%' OR DIAG LIKE 'DI110%' OR DIAG LIKE 'DI130%' OR DIAG LIKE 'DI132%' OR DIAG LIKE 'DI255%' OR DIAG LIKE 'DI420%' OR DIAG LIKE 'DI425%' OR DIAG LIKE 'DI426%' OR DIAG LIKE 'DI427%' OR DIAG LIKE 'DI428%' OR DIAG LIKE 'DI429%' OR D

### Alder/korrekt århundrede udfra cpr nummer

De fleste, der arbejder med registre eller databaser, står ofte med problemstillingen, at alder er uoplyst, medens cpr-nummer er kendt. Hvordan regner man den ud? Følgende regel er gældende: Hvis syvende ciffer er 0, 1, 2 eller 3 er man født i det 20. århunderede (1900-tallet) Ligeledes, hvis syvende ciffer er 4 eller 9, og årstallet (femte og sjette ciffer) er større end eller lig 37. Endelig er man født i det 19. århundrede (1800-tallet) hvis syvende ciffer er 5, 6, 7 eller 8 og årstallet er større end eller lig 58. Nedenfor finder du eksempel i SAS kode: En lille makro, der udover fødselsdato også udregner køn samt den præcise alder givet datovariabel. Kilde: Opbygning af CPR nummeret, cpr.dk proc format library=work; value gender 0="Female" 1="Male" ; run; %macro agefromCPR(cpr,datevar=inddto,birthvar=birth,agevar=age); dy_temp=input(substrn(&cpr,1,2),2.); mt_temp=input(substrn(&cpr,3,2),2.); yr_temp=input(substrn(&cpr,5,2),