Minimal effort method for integrating C++ functions into Mathematica

Question

Update: While at the time of writing the question loading DLLs with .NET/Link seemed easier, now I always use LibraryLink, which I recommend to anyone with a similar problem!

---

As of Mathematica 8, what is the minimal effort way to integrate an existing C++ function into Mathematica?

I think we have these:

• MathLink (it was quite long ago I used it last time)
• communication through pipes/files (Import is slow, ReadList not so much)
• LibraryLink (??)
• Mathematica 8's C code generation features (??) _{(apparently this is not relevant)}

The keyword here is minimal effort. Which is the most convenient way, including the learning curve for the particular method?

_{I'm mainly interested in doing it on Windows. My particular function computes a long list of real numbers. The function input consists only of a few scalars (int & double).}

_{An answer of the type "You'll likely spend the least time if you use technology X" is useful. A concrete example of how to do it using a small function is even better.}

---

Answers

Each and every answer I got is great, and each demonstrates a different technology. It's impossible to choose a definitive one. Here's a small "table of contents" for them:

• MathLink
• Read result from pipes, pass arguments on command line to separate process
• LibraryLink
• Using .NETLink for loading functions from arbitrary DLLs

A presentation about the topic from the Wolfram Technology Conference is here (CDF file).

---

Experiences

I used two of these methods: the .NETLink method and Library Link. These are my personal experiences:

• .NETLink This was easy to set up, however, I'd recommend it only for when you already have a compiled DLL from which you need to load functions (e.g. do something like this). If you compile your own DLL: once the DLL is loaded, it is locked and cannot be overwritten. A quick and dirty way to allow overwriting it is realoding NETLink with ReinstallNET[].

  Advantages: Very quick and easy, only Mathematica code is needed if the functions are already compiled into a DLL.

  Disadvantages: Windows-only, does not parallelize from within Mathematica, and the calculation will not be interruptible.

• Library Link It's much easier to use this than what it looks like at first. Less setup is needed than in the case of MathLink, and compilation is automated from within Mathematica.

  Advantages: Also quite easy, but both Mathematica and C code are needed. It is parallelizable, and easily made interruptible. The library can be unloaded (for recompilation) using LibraryUnload.

Answer 1

Update: I posted a tutorial series on LibraryLink at Wolfram Community.

Update: LibraryLink example code is also available on this GitHub repository.

---

Here is my pitch to use LibraryLink, which is a really nice new technology in version 8. I am not going to pretend this is easy by any stretch of the imagination, because it involves a decent amount of knowledge of both Mathematica and C compilers. In this particular case I am using Visual Studio C++ Express plus the Microsoft Windows SDK 7.1. For the record, I had quite a bit of help from Joel Klein with this answer.

LibraryLink is set up to find most compilers on your system, but in my case I had to restart after installing the above tools (although looking back I think restarting my open frontend might have done the trick as well).

There are several examples of how to use LibraryLink in the documentation, but this example I wrote from scratch. All LibraryLink C code is compiled with the CreateLibrary function (which is located in the CCompilerDriver package.

Mathematica side

I'll skip the source for now, and focus on the Mathematica commands. First we loads the package to run the library link utility functions:

Needs["CCompilerDriver`"]

Next, we load the source file and create a library (.dll) from it:

myLibrary = 
  CreateLibrary[{"c:\\users\\arnoudb\\myLibrary.c"}, "myLibrary", "Debug" -> False];

Next, we load the one function defined in this library:

myFunction = LibraryFunctionLoad[myLibrary, "myFunction", {{Real, 2}}, {Real, 2}];

Next we use this function:

myFunction[{{1, 2, 3, 4}, {5, 6, 7, 8}}]

which returns the square of every matrix entry: {{1., 4., 9., 16.}, {25., 36., 49., 64.}}.

Finally, we unload the library (needed if we make changes to the source file and reload everything above):

LibraryUnload[myLibrary];

C side

There is a lot of boiler plate magic, to make things 'work'. These first four line need to always be included:

#include "WolframLibrary.h"
DLLEXPORT mint WolframLibrary_getVersion(){return WolframLibraryVersion;}
DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {return 0;}
DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {}

This is the actual function that you write. The function header is always the same, you get the actual function arguments from the MArgument_* api functions:

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

 int err; // error code

 MTensor m1; // input tensor
 MTensor m2; // output tensor

 mint const* dims; // dimensions of the tensor

 mreal *data1; // actual data of the input tensor
 mreal *data2; // data for the output tensor

 mint i; // bean counters
 mint j;

This gets the input tensor and dimensions and sets up the output tensor and actual data:

 m1 = MArgument_getMTensor(Args[0]);
 dims = libData->MTensor_getDimensions(m1);
 err = libData->MTensor_new(MType_Real, 2, dims,&m2);
 data1 = libData->MTensor_getRealData(m1);
 data2 = libData->MTensor_getRealData(m2);

The actual interesting stuff, this squares each element:

 for(i = 0; i < dims[0]; i++) {
  for(j = 0; j < dims[1]; j++) {
   data2[i*dims[1]+j] = data1[i*dims[1]+j]*data1[i*dims[1]+j];
  }
 }

This set the return value (and yes, you don't want to return the 'err' value here):

 MArgument_setMTensor(Res, m2);
 return LIBRARY_NO_ERROR;
}

Answer 2

On Windows, C/C++ functions that have been compiled into DLLs can be accessed reasonably easily using NETLink. Let's say we have the following C++ DLL definition:

#include "stdafx.h"

BOOL APIENTRY DllMain(HMODULE module, DWORD reason, LPVOID reserved) {
    return TRUE;
}

extern "C" __declspec(dllexport)
void helloMma(double a, double b, int n, double m[]) {
    for (int i = 0; i < n ; ++i) {
        m[i] = a * i + b;
    }
}

We can import this function into Mathematica as follows:

Needs["NETLink`"]

$dllPath = "C:\\some\\path\\to\\hellomma.dll";

helloMma = DefineDLLFunction[
  "helloMma", $dllPath, "void", {"double", "double", "int", "double[]"}
]

Alas the calling sequence for this function is complicated slightly by the fact that it returns its result by destructively overwriting a statically allocated array of doubles (not an uncommon occurrence). To do this from Mathematica, we must pre-allocate the result vector using NETNew:

In[23]:= NETBlock@Module[{n, result}
         , n = 10
         ; result = NETNew["System.Double[]", n]
         ; helloMma[3, 5, n, result]
         ; NETObjectToExpression[result]
         ]
Out[23]= {5., 8., 11., 14., 17., 20., 23., 26., 29., 32.}

Note that all the usual Bad Things would happen if the pre-allocated buffer were overrun. NETBlock is used to ensure that the storage allocated for the buffer is released when we are done with it.

I will point out a couple of "gotchas". Make sure that the DLL is compiled as 32-bit or 64-bit to match the version of Mathematica that you are running. Also, note that the exported DLL function in the example is declared as extern "C". This prevents the C++ compiler from "mangling" the name and makes it easier to reference in the DefineDLLFunction declaration (in this case the Visual Studio C++ mangled name was ?helloMma@@YAXNNHQEAN@Z).

Answer 3

Presuming that your c++ code is already written, then I don't know how the code generation feature would be helpful. That said, in order of simplicity, I would have Get, ReadList, Import, and both LibraryLink and MathLink.

---

Get and ReadList are by far the simplest. Provided that your c++ program outputs to stdout (std::cout), then it is simply

val = (<<"!command")

or,

ReadList["!command", Number (*or some other specifier*)]

If you need to control your external program more directly, or pass additional data to it (like run your function multiple times), then this method is more difficult and may require you to open pipes to communicate through (1, 2).

---

Import would require you to conform to some data format, which while very doable, is not as simple as just putting numbers onto stdout.

---

I have never used LibraryLink, but a quick perusal of the docs implies that it is a simplification of the MathLink interface. More specifically, it is easier to create a function that can communicate with Mathematica, but the intricacies of sending data back and forth remain. To return a List, though, is not to bad, and may well be worth it. The catch, of course, is that you need to provide a c interface to your function, which differs from the MathLink form, but is still necessary. So, I'm not sure there is an advantage here for one or the other.

Answer 4

I suggest to use MathLink, which you can automate using the CCompilerDriver`. This is a safe alternative, since you won't crash the kernel if your code crashes. Once tested, this should not be hard to convert to library link. As an explicit example, consider a function which receives a list of integers and squares it. First, here is a function to create the boilerplate code:

makeMLinkCodeF = 
  StringJoin[
     "#include \"mathlink.h\"", "\n", ##, "\n",
     "
     #if WINDOWS_MATHLINK

     #if __BORLANDC__
     #pragma argsused
     #endif

     int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE 
           \ hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
     {
     \tchar  buff[512];
     \tchar FAR * buff_start = buff;
     \tchar FAR * argv[32];
     \tchar FAR * FAR * argv_end = argv + 32;

     \thinstPrevious = hinstPrevious; /* suppress warning */

     \tif( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;
     \tMLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
     \treturn MLMain( (int)(argv_end - argv), argv);
     }

     #else

     int main(int argc, char* argv[])
     {
     \treturn MLMain(argc, argv);
     }

     #endif"] &;

It is pretty ugly, but you only need to define it once. Perhaps, a better way would be to rewrite it using Symbolic C, but I had no time to do it. Now, here is the code for our function:

code = 
"
 extern void squareList(int * data, long len);

 void squareList(int * data, long len){
    int *d = data;
    MLPutFunction(stdlink,\"List\",len);
    while(d-data<len){
        MLPutInteger(stdlink,(*d) * (*d));
        d++;
    }
 }
";

And here is the template:

template = 
StringReplace[#,"\n"~~Whitespace:>"\n"]&@
"
void squareList P((int *, long));

:Begin:
:Function:       squareList
:Pattern:        squareList[data_List]
:Arguments:      { data }
:ArgumentTypes:  { IntegerList }
:ReturnType:     Manual
:End:

";

(StringReplace was needed since otherwise the template would not copy-paste correctly from SO - it would contain spaces on the left, which would prevent us from creating an executable. You don't have to do this in the interactive session).

We must now load the compiler driver:

Needs["CCompilerDriver`"];

The following code prepares our stuff:

fullCCode = makeMLinkCodeF[code];
projectDir = "C:\\Temp\\MLProject";
If[! FileExistsQ[projectDir], CreateDirectory[projectDir]]
pname = "squareList";
files = 
 MapThread[
    Export[FileNameJoin[{projectDir, pname <> #2}], #1, "String"] &, 
    {{fullCCode, template}, {".c", ".tm"}}];

Now we attempt to create an executable:

In[19]:= exe=CreateExecutable[files,pname]
Out[19]= C:\Users\Archie\AppData\Roaming\Mathematica\SystemFiles\LibraryResources
  \Windows-x86-64\squareList.exe

We now install it:

In[20]:= Install[exe]
Out[20]= LinkObject["C:\Users\Archie\AppData\Roaming\Mathematica\SystemFiles\
LibraryResources\Windows-x86-64\squareList.exe",10,8]

And finally use it:

In[21]:= squareList[Range[5]]
Out[21]= {1,4,9,16,25}


In[23]:= Uninstall[exe]
Out[23]= C:\Users\Archie\AppData\Roaming\Mathematica\SystemFiles\LibraryResources
\Windows-x86-64\squareList.exe

There were a number of steps here, but all of them are within Mathematica, and most of them are common for all functions and can be factored away. If you end up doing lots of coding, I'd create a Mathematica-based DSL which would generate Symbolic C for your functions and also your templates. If / when time permits, I intend to improve this answer by automating the procedure more, and by replacing the fragile string-based pieces with functions to generate them.

Answer 5

In terms of loading C++ functions, I would strongly suggest LibraryLink. It is a great tool except that it requires you to write sometimes intimidating C code.

To make LibraryLink easier to use, I have developed a package called wll-interface, available in this github repository. It is a header-only library written in C++, doing only one thing—automatically generate the code that interact with LibraryLink.

---

A simple example

Here I give an example of how wll-interface simplifies the process of loading a C++ function; this example can also be found here.

Let's say that we have a C++ function multiply and we want to load it into Mathematica.

double multiply(double x, int y)
{
    return x * y;
}

Instead of calling LibraryLink APIs, all we need to do with the help of wll-interface is 1) include the header file, and 2) use DEFINE_WLL_FUNCTION macro. This is what the source code would look like:

src = "
#include \"wll_interface.h\"  // include the header file

double multiply(double x, int y)
{
    return x * y;
}

DEFINE_WLL_FUNCTION(multiply)  // defines wll_multiply
";

Then we create a shared library from the code. You need to replace <path-to-wll_interface.h> below with the directory that contains the header file wll_interface.h so that the compiler can find it.

Needs["CCompilerDriver`"];
mylib = CreateLibrary[src, "wll_multiply"(* name of the library, can be anything *), 
  "IncludeDirectories" -> {"<path-to-wll_interface.h>"}, 
  Language -> "C++", "CompileOptions" -> "-std=c++17"]

Finally we load the function from the shared library we just created. Note that we are loading wll_multiply—a function that is automatically generated.

multiply = LibraryFunctionLoad[mylib, "wll_multiply", {Real, Integer}, Real];

Now multiply can be used just like a normal function!

multiply[2.33, 5]

---

Manipulating arrays (MTensor)

The template class wll::tensor<type, rank> encapsulates LibraryLink MTensor, making operations on arrays much easier. In this example, I will show you how to pass a reference to an array through LibraryLink and sort the array.

wll::tensor has member functions begin and end to get the iterators, similar to STL containers. Declaring the type of x to be wll::tensor<double, 1>& means that this list of doubles will be passed by reference and no copy will be performed.

src = "
#include <algorithm>
#include \"wll_interface.h\"

void sort(wll::tensor<double, 1>& x)
{
    std::sort(x.begin(), x.end());
}
DEFINE_WLL_FUNCTION(sort)
";

Now we compile the library as before and load the library function wll_sort

Needs["CCompilerDriver`"];
mylib = CreateLibrary[src, "wll_sort",
  "IncludeDirectories" -> {"<path-to-wll_interface.h>"}, 
  Language -> "C++", "CompileOptions" -> "-std=c++17"]
sort = LibraryFunctionLoad[mylib, "wll_sort", {{Real, 1, "Shared"}}, "Void"];

Note that the type of argument is specified as {Real, 1, "Shared"}, where "Shared" means pass by reference here.

This library function works by directly modifies values of the variable in the Wolfram Kernel:

a = RandomReal[1., 10];
sort[a]   (* sort the array; there is no return value *)
a         (* the array has been sorted! *)

---

Sparse arrays (MSparseArray)

In this example, we will use wll-interface to help us generate identity matrices as sparse arrays. The template class wll::sparse_array<type, rank> corresponds to MSparseArray in LibraryLink. Its member function operator() can be used to extract the elements (wll::tensor also have this function).

src = "
#include \"wll_interface.h\"

wll::sparse_array<double, 2> sparse_identity(size_t size)
{
    wll::sparse_array<double, 2> identity({size, size});
    for (size_t i = 0; i < size; ++i)
        identity(i, i) = 1.0;
    return identity;
}
DEFINE_WLL_FUNCTION(sparse_identity)
";

Here identity(i, i) points to an element on the main diagonal. Following the convention of C/C++, indices in wll-interface are zero-based.

Compile and load the function:

Needs["CCompilerDriver`"];
mylib = CreateLibrary[src, "wll_sparse_identity",
  "IncludeDirectories" -> {"<path-to-wll_interface.h>"}, 
  Language -> "C++", "CompileOptions" -> "-std=c++17"]
identity = LibraryFunctionLoad[mylib, "wll_sparse_identity", {Integer}, LibraryDataType[SparseArray]];

identity[10]  (* gives a 10x10 idensity matrix *)

Answer 6

There is a package called LTemplate that automates writing some of the boilerplate code for LibraryLink:

• How to simplify writing LibraryLink code?

I consider this less effort than writing standard LibraryLink code. In this sense it is a fitting answer for this question.

However, I do recommend familiarizing yourself with the standard way of using LibraryLink before trying LTemplate.

Disclosure: I am the author of LTemplate.

Answer 7

You may want to look at this presentation: Integrating C and Mathematica.

In the past, I have found using .NET/Link to be the easiest. You can call C DLL's very easily on Windows, without the need for templates as in MathLink.