Linking FORTRAN with Librarylink

Question

FORTRAN code can be called using MathLink or .NET/Link (see the link for a worked examples).

But as mentioned in a talk by T.Gayley and J.Klein in Wolfram Technology Conference 2011, LibraryLink, which is new in version 8, is much faster. In the presentation, they gave a comparison of these three different approaches. (see also here a comparison by @Szabolcs). I just copy the pros and cons of LibraryLink from that presentation:

---

• Pros

  • Extremely fast, comparable to internal implementation in Mathematica itself.

  • Memory for large arrays can be shared with the kernel.

  • Code generation tools in Mathematica can be used to simplify creating and compiling libraries.

  • Probably a better alternative than installable MathLink programs for most users who want to call C functions.

• Cons

  • Need to write specialized C code for every function you want to call.

  • Without special (and sometimes inconvenient) effort by the author, the kernel is not interruptible while calling a library function. This means that functions cannot be aborted while in progress, and the front end will not be responsive to most buttons, the help system, etc.

---

The advantages arouse my interest. But I have little experience with C programming... :-( So I want to know how can one call fortran functions and subroutines using Librarylink? Is there any worked example?

Practically, I want to use it to call a cernlib package MINUIT written in FORTRAN 77, see my question here. MINUIT is a complicated package for function minimization and error analysis, and it has been widely used for may years, especially in the fields of nuclear and particle physics. Here is a link to its reference manual.

Answer 1

Here are three very simple examples to show how to call a Fortran subroutine using LibraryLink. First the subroutine is compiled into object file. Then a wrapper is used to call the Fortran subroutine and compiled into dynamic library. At the end, the library is loaded into Mathematica and run. In the examples Mathematica Version 8 is used.

---

FIRST EXAMPLE

add two integers

Fortran subroutine

!fadd.f90
subroutine add(a,b,sum)
    implicit none
    integer a,b,sum
    sum=a+b
    return
end subroutine

C Wrapper

//MMA.cc
//Link directly to fortran object file
#include "WolframLibrary.h"

DLLEXPORT mint WolframLibrary_getVersion(){
    return WolframLibraryVersion;}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
    return 0;}
extern "C" {
    void add_(mint* a, mint* b, mint* sum);} //declare fortran subroutine

EXTERN_C DLLEXPORT int add(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){
    mint I0;
    mint I1;
    mint sum;
    I0=MArgument_getInteger(Args[0]);
    I1=MArgument_getInteger(Args[1]);

    add_(&I0,&I1,&sum);//call fortran subroutine

    MArgument_setInteger(Res,sum);
    return LIBRARY_NO_ERROR;
}

Compile and Run

First compile fortran subroutine:

gfortran -c fadd.f90

Then create dynamic library using the CCompilerDriver package:

Needs["CCompilerDriver`"];
CreateLibrary[{"MMA.cc", "fadd.o"}, "myadd", "Debug" -> True, "TargetDirectory" -> "."]

Load and run:

add = LibraryFunctionLoad["./myadd", "add", {Integer, Integer}, Integer]
add[2, 2]
(*LibraryFunction[<>,add,{Integer,Integer},Integer]*)
(*4*)

---

SECOND EXAMPLE

add two vectors (code modified from similar question at here).

Fortran subroutine

!addvec.f90

subroutine addVec(a,b,n)
implicit none
integer n
real(8),dimension(n)::a,b

a(:)=a(:)+b(:)

return
end subroutine addVec

C Wrapper

//MMA.cc

#include "WolframLibrary.h"
#include "WolframCompileLibrary.h" 

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
  return 0;
}
extern "C" {
  void addvec_(mreal a[], mreal b[], mint* n);
} 


EXTERN_C DLLEXPORT int addvec(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

  MTensor ta;
  MTensor tb;

  ta=MArgument_getMTensor(Args[0]);
  tb=MArgument_getMTensor(Args[1]); 
  addvec_(MTensor_getRealDataMacro(ta),MTensor_getRealDataMacro(tb),MTensor_getDimensionsMacro(ta));

  MArgument_setMTensor(Res,ta);
  return LIBRARY_NO_ERROR;
}

Compile and Run

Needs["CCompilerDriver`"]
CreateLibrary[{"MMA.cc", "addvec.o"}, "myadd", "Debug" -> True, "TargetDirectory" -> "."]
addVec = LibraryFunctionLoad["./myadd", "addvec", {{Real, 1}, {Real, 1}}, {Real, 1}]

addVec[{1.4, 2.7, 3.9}, {5.2, 6.7, 7.1}]
addVec[{1.1, 2.2, 3.2, 2.2}, {4.4, 2.2, 3.3, 5.5}]
(*{6.6, 9.4, 11.}*)
(*{5.5, 4.4, 6.5, 7.7}*)

---

Third EXAMPLE

invoke Lapack to inverse a general matrix (code modified from here).

C Wrapper

//MMA.cc
#include "WolframLibrary.h"
#include "WolframCompileLibrary.h"

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
  return 0;
}

extern "C" {
  // LU decomoposition of a general matrix
  void dgetrf_(mint* M, mint *N, double* A, mint* lda, mint* IPIV, mint* INFO);    
  // generate inverse of a matrix given its LU decomposition
  void dgetri_(mint* N, double* A, mint* lda, mint* IPIV, double* WORK, mint* lwork, mint* INFO);
}

EXTERN_C DLLEXPORT int lpkInverse(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

  MTensor ta=MArgument_getMTensor(Args[0]);
  mint N=*MTensor_getDimensionsMacro(ta);    
  double* A=MTensor_getRealDataMacro(ta);

  mint IPIV[N+1];
  mint LWORK = N*N;
  double WORK[LWORK];
  mint INFO;

  dgetrf_(&N,&N,A,&N,IPIV,&INFO);
  dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);

  MArgument_setMTensor(Res,ta);    
  return LIBRARY_NO_ERROR;    
}

Compile and Run

Needs["CCompilerDriver`"]
$CCompiler = {"Name" -> "g++", "Compiler" -> CCompilerDriver`GenericCCompiler`GenericCCompiler, "CompilerInstallation" -> "/usr/bin", "CompilerName" -> "g++"};
CreateLibrary[{"MMA.cc"}, "myadd", "Debug" -> True, "TargetDirectory" -> ".", "CompileOptions" -> "-llapack"]
inverse = LibraryFunctionLoad["./myadd", "lpkInverse", {{Real, 2}}, {Real, 2}]


inverse[{{1, 2}, {3, 4}}]
Inverse[{{1, 2}, {3, 4}}] // N
(*{{-2., 1.}, {1.5, -0.5}}*)
(*{{-2., 1.}, {1.5, -0.5}}*)


a = Table[Table[RandomReal[{0., 10.}], {10}, {10}], {100}];
bM = Inverse /@ a; // AbsoluteTiming
bl = inverse /@ a; // AbsoluteTiming
(*{0.005429, Null}*)
(*{0.001198, Null}*)

bM - bl // Abs // Max
(*1.30562*10^-13*)

I'm still in learning LibraryLink so if there are any mistakes please don't hesitate to point them out. Hope it will help.