Passing Mathematica function via MathLink or LibraryLink to c (or c++)

Question

I want to numerically solve the heat equation. To do so I wrote a program in C++ - It is a BVP $(\alpha T +\kappa \frac{\partial T}{\partial r}=f(t))$, where $f(t)$ is an arbitrary function, whose address I pass to the constructor of the struct BoundaryCondition

struct BoundaryCondition
{
  enum class Side {LEFT,RIGHT};
  double alpha; 
  double kappa;
  double (*f) (double);
  Side side;
  BoundaryCondition(Side,double,double,double (*)(double));
};

Now I rewrite my code to use Mathematica Interface.

My plan is to write the functions for the left and right boundary conditions in Mathematica, pass these functions(by name or by address???) and the initial data(array of real numbers) to my C++ code. Then I will just call these functions in an iterative algorithm. And finally - flush the results back to Mathematica for post-processing (graphics, another calculation, etc.)

Do you have any idea how to pass an arbitrary (unknown in advance) function (not just the usual data) to my, written in C++, library as an argument for further usage in the C++ code ?

Answer 1

This is too long for a comment, but let me share some ideas. For the second part of the answer, I'm very skeptical whether it will work, but if you are investigating into the issue, you probably want to check this path.

First of all:

where f(t) is an arbitrary function

For this general case you won't be able to always create a compiled function, because only a small part of Mathematica constructs can be compiled down into a library which doesn't need a mathKernel for evaluation. Please see this thread for more information about what is compilable in Mathematica.

In the case your function is really arbitrary, then I see only one way and this is to use MathLink to communicate with the Kernel from within you C++ function. With this you can send the values to Mathematica, let it evaluate the function f and read back the result.

To say it plain: this is not recommended in your situation, because as you try to solve a pde numerically, you right hand side is going to be evaluated very often. I'm afraid the overhead of the communication will destroy any speed advantage.

Let's assume f[x] is friendly enough so that you could create a compiled function, which itself doesn't run evaluations through the kernel by MainEvaluate calls. Then, there might be an idea which you could try:

Before you call your C++ library pde-solver, you compile f[x] to "C" which creates a shared library. You have to make this in a way that you know how the library and function name. These information can then be send to your C++ library (or they are fixed). Within your C++ library, you dynamically load the created dll and retrieve the function pointer and after that, you might be able to call your function f.

To give you a start, first I would look how the library function created by Mathematica looks. For this purpose the CCodeGenerator` package is required:

<< CCodeGenerator`
f = Compile[{{x, _Real, 0}}, x^2];
CCodeStringGenerate[f, "f"]

I will not post the whole output, but you see that there is a Initialize_f function, which is important to set up the global funStructCompile. You surely want to call this before doing anything else.

Your function f in C-code looks as follows:

DLLEXPORT int f(WolframLibraryData libData, mreal A1, mreal *Res)
{
  mreal R0_ 0;
  mreal R0_ 1;
  R0_ 0 = A1;
  R0_ 1 = R0_ 0 * R0_ 0;
  *Res = R0_ 1;
  funStructCompile->WolframLibraryData_cleanUp(libData, 1);
  return 0;
}

You see, you only need to provide the argument as mreal and a pointer for the result.

For the creation of the dll, you don't need to inspect the C-code every time. I just showed it so that you see how the naming happens inside the library. What you need is

LibraryGenerate[f, "f"]

which returns the place where the created dll is.

I'm not sure whether I mentioned everything important. If something is unclear, please ask, because I would really like to see whether this works.

One helpful link: Dynamically load a function from a DLL

Answer 2

Since Mathematica 10.0, you can use library callback functions in LibraryLink, which are designed precisely for this purpose.