Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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
1  
@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 7:33

1 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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.