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.

Keep in mind that one possible disadvantage of LibraryLink is that if your code crashes the kernel will crash, and you may then lose your results. I would still first use MathLink to prototype the solution, and then port to LibraryLink for speed once you are sure that your code does not crash often. — Leonid Shifrin, May 18 '12 at 13:44
@Leonoid, I tried Mathlink first. For some unknown reason, it didn't produced a link... If would be nice if you could have a look at my quesion here. — unstable, May 18 '12 at 13:48
It would be nice if we could see the LibraryLink version here. I'm just wondering, would this execute faster than the MathLink version of the program? — jmlopez, May 18 '12 at 17:30
@jmlopez the NETLink version seems to be faster than Mathlink. I'll have a second look at the mathlink problem following your answer. My NETLink version of code finishes a realistic three-parameter fit to 20 data points in less than 0.4 second (no integral). The time of install a mathlink is typically more than 20 seconds in my experience. — unstable, May 18 '12 at 18:50
Does anyone happen to know if the LibraryLink API is a superset of the MathLink API insofar as accessing frontend internals is concerned? I'm looking for anything that might expose hooks that can be set much like in the emacs model for editor customization. — StackExchanger, May 20 '12 at 07:33
@unstable, were you able to get the Mathlink version working? If so, have you started looking into the LibraryLink version? — jmlopez, May 22 '12 at 02:18
@jmlopez, since now NETLink is working rather well and quickly, I'm busy writing my fortran code and calculating -- into late night almost everyday... and have not tried the mathlink yet, neither librarylink. I'll do that later. — unstable, May 22 '12 at 08:40
@StackExchanger Neither MathLink, not LibraryLink provides any access to Front End internals. LibraryLink gives a little bit of access to kernel internals: direct access to in-memory numerical arrays. — Szabolcs, Feb 14 '13 at 17:59
@unstable Are you sure NETLink is faster than MathLink? That's surprising to me as NETLink also uses MathLink in the background. (But of course there's more than one way to use MathLink) — Szabolcs, Feb 14 '13 at 18:02
@unstable Have you already solved this problem? I see you have a StackOverflow answer which explains how to do this using MathLink. The procedure should be similar with LibraryLink: you will need a little bit of C code which interfaces with LibraryLink and the FORTRAN library. Fortunately C is flexible and supports various calling conventions, so FORTRAN routined can be called form C directly (with care), just as you showed in your StackOverflow answer I linked to. The procedure here is the same: write a C function which work with LibraryLink and call another ... — Szabolcs, Feb 14 '13 at 18:05
... FORTRAN function the appropriate way. (I understand that you're not familiar with C so this would be some trouble for you to figure out, I'm just trying to give some guidance.) — Szabolcs, Feb 14 '13 at 18:06
@Szabolcs, not worked on that so far. Now I use NETLink, and it's fast. Not sure whether or not faster than MathLink. I remember that it always took me some seconds to establish a link in Mathmetica 8 in my old laptop. With my new laptop and Mathematica 9, the connection time seems order-of-magnitude shorter... — unstable, Feb 19 '13 at 09:28

score 38 · Accepted Answer · edited May 23 '17 at 12:35

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.

thanks for a nice tutorial; FYI the second example works for me but not the first one. — chris, Oct 18 '15 at 18:16
Hi, great answer. I am new to librarylink, really want to master it. Right now, I just want to call lapack in mma to speed up inverse in my program. But I am in windows and using intel fortran compiler ifort. Is it possible to call intel MKL library using library link? — matheorem, Dec 22 '15 at 16:05

Linking FORTRAN with Librarylink

1 Answers1

Linked