45
$\begingroup$

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.

Improve this question
$\endgroup$
12
  • 2
    $\begingroup$ 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. $\endgroup$ Commented May 18, 2012 at 13:44
  • 1
    $\begingroup$ @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. $\endgroup$ Commented May 18, 2012 at 13:48
  • $\begingroup$ 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? $\endgroup$ Commented May 18, 2012 at 17:30
  • $\begingroup$ @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. $\endgroup$ Commented May 18, 2012 at 18:50
  • $\begingroup$ 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. $\endgroup$ Commented May 20, 2012 at 7:33

1 Answer 1

Reset to default
39
+500
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ thanks for a nice tutorial; FYI the second example works for me but not the first one. $\endgroup$ Commented Oct 18, 2015 at 18:16
  • $\begingroup$ 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? $\endgroup$ Commented Dec 22, 2015 at 16:05

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.