How to call the C++ by the new LibraryFunctionDeclaration?

Question

When we arrive at V13.1, we have a great LibraryFunctionDeclaration now, which can call C function in MMA as this article:

Needs["CCompilerDriver`"]
lib = CCompilerDriver`CreateLibrary["#include \"WolframLibrary.h\"
#include \"WolframCompileLibrary.h\"

template<typename T>
void selection_sort(T* a, int low, int high)
{
int min;
    for (int i = low; i <= high; i++)
    {
    min = i;
    for (int j = i + 1; j <= high; j++)
    {
    if (a[j] < a[min])
    min = j;
    }
    if (min != i)
    {
    auto temp = a[min];
    a[min] = a[i];
    a[i] = temp;
    }
    }
}

EXTERN_C DLLEXPORT MTensor _selectionSort_Cpp(MTensor in)
{
auto in_data = (mreal*)in->data;
(in->refcount)++;
selection_sort(in_data, 0, in->nelems - 1);
return in;
}", "selectionSort", Language -> "C++"]
dec = LibraryFunctionDeclaration[
   selectionSort -> "_selectionSort_Cpp", 
   lib, {"PackedArray"["Real64", 1]} -> "PackedArray"["Real64", 1]];
fc = FunctionCompile[dec, 
  Function[Typed[n, "PackedArray"["Real64", 1]], selectionSort[n]]]

Then we can use C function(selection_sort) directly in MMA now:

lis = RandomReal[{-10, 10}, 100000];
fc[lis]; // Timing

{1.54688, Null}

It is awesome. But since we can use the Language -> "C++". Why not use C++ grammar? Such as we have a function sort_add2 like in the following:

#include<iostream>
#include<vector>
#include <algorithm>

using namespace std;

vector<int> sort_add2(vector<int> vec) {
    sort(vec.begin(), vec.end());
    for (int &i : vec) {
        i=i+2;
    }
    return vec;
}

int main()
{
    vector<int> vec = {5,9,20,3};
    vector<int> new_vec = sort_add2(vec);
    for(int i : new_vec){
        cout << i << endl;
    }

    return 0;
}

How to call that sort_add2 directly in MMA 13.1? Is it enough to construct a structure with MMA's existing compiled type so that MMA can recognize the input and output of sort_add2? Is possible?

Answer 1

This post is written with Mathematica version 13.1 and is subject to change in the future.

Mathematica's C Types does not include vector type but includes CArray which can be converted to C vector type easily.

Compiling, in general, and in this case, will change symbol names that's why you used EXTERN_C DLLEXPORT which are extern "C" and __declspec(dllexport) respectively to make sure they're exported and their names won't change. So I added a little bit of interface to make the code work.

1. Let's assume our code looks like this:

cppCode = "
#include <vector>
#include <algorithm>

using namespace std;

// accept a reference to change the original data
void sort_add2(vector<int>& vec) {
    
    sort(vec.begin(), vec.end());
    for (int &i : vec) {
        i=i+2;
    }
}

extern \"C\" {
    __declspec(dllexport) int* sort_add2_interface(int* arr,int len){
        
        // build a vector from our array
        vector<int> temp(arr, arr + len);

        // do the operation
        sort_add2(temp);

        // return the first element address
        return &(temp[0]);
    }
}
";

2. Now build the code:

lib = CCompilerDriver`CreateLibrary[cppCode, "selectionSort", 
  Language -> "C++"]

3. Declare the function:

declaration = 
 LibraryFunctionDeclaration["sort_add2_interface", 
  lib, {"CArray"::["CInt"], "CInt"} -> "CArray"::["CInt"]]

4. Load the exported library:

LibraryLoad[lib];

5. On the Mathematica side, we convert our list to CArray using CreateTypeInstance, then call the library and build the result from reading the output with FromRawPointer (it get casts to Integer64).

fn = FunctionCompile[declaration, 
  Function[Typed[arr, "ListVector"::["Integer32"]], 
   Module[{temp, len, result},

    len = Cast[Length[arr], "CInt", "BitCast"];

    temp = 
     LibraryFunction["sort_add2_interface"][
      CreateTypeInstance["Managed"::["CArray"::["CInt"]], arr], len];

    Table[Cast[FromRawPointer[temp, i - 1], "Integer64"], {i, len}]

    ]
   ]
  ]

Testing:

fn[{2, 3, 1}]

(* Out: {3, 4, 5} *)

Remember that if you have a library that contains sort_add2 and is exposed, you can build another library that acts as an interface. In this case, I combine the two but they can be separated.

If they are separated, you have to include their declaration to your code and use CCompilerDriver`CreateLibrary's "Libraries" option to make it work.

For example, if you have a lib1.dll (with a lib1.lib) which contains sort_add2, you code will look like this:

cppCode2 = "
#include <vector>
#include <algorithm>

using namespace std;

// declaration
extern \"C\" void sort_add2(vector<int>& vec);

extern \"C\" __declspec(dllexport) int* sort_add2_interface(int* \
arr,int len){
    vector<int> temp(arr,arr+len);
    sort_add2(temp);
    return &(temp[0]);
}
";

You would use this code to compile:

lib = CCompilerDriver`CreateLibrary[cppCode2, "library_interface", 
  Language -> "C++", "Libraries" -> {"C:\\PATH_TO_LIBRARY\\lib1.lib"}]

And the rest is similar (follow from step 3, in step 4 you need to load both libraries).

---

I was more curious about the performance, so here is the result (assuming data = RandomInteger[{0, 2^31 - 1}, 10000]):

Code	Timing	Memory Usage
Sort[data]+2	0.00025	160,456
fn[data]	0.0013 (5x slower)	281,160

I decide to go deeper and this is the result of different parts of calling fn1:

Code	Timing (~)
building vector (CreateTypeInstance)	0.00083
calling sort_add2	0.00048
building the result (Table[ ... FromRawPointer, ... ])	0.00001

The newly added feature has much more power and along with the LibraryLink which can give direct access to the data (and I assume is faster), Mathematica become more connected than ever.

There are many things that I don't and I'm sure the code above can improve further. (I'll welcome any suggestion/correction)

Update 1

Thanks to @yode's comment and this post using "PackedArray" instead of "ListVector" will boost the speed by a factor of 2, so the new code will look like this:

fn2 = FunctionCompile[declaration, 
 Function[Typed[x, "PackedArray"["MachineInteger", 1]],
  Module[{in, len, out},
   len = Cast[Length[x], "CInt", "BitCast"];
   in = CreateTypeInstance["CArray"["CInt"], len];
   Do[ToRawPointer[in, i - 1, Cast[x[[i]], "CInt", "BitCast"]], {i, 
     len}];
   out = LibraryFunction["sort_add2_interface"][in, len];
   Table[Cast[FromRawPointer[out, i - 1], "Integer64"], {i, len}]]]]

In pursuing performance I did try sort_add2 with the i64 type, so we might gain speed by removing all the Cast but it didn't have an effect.

My next attempt was building a library and hoping by using LibraryFunctionLoad, we could gain speed which didn't happen (both Automatic and "Shared" memory management types). It was around ~20% slower than fn2. Interestingly whether it was "Shared" (in-place operation) or it was a copy of the original data, the speed didn't have a difference.

Answer 2

This is my program, in fact most of it comes from an unknown author. And it exactly is what I'm after. For more details you need to see follow this article for updates

---

Needs["CCompilerDriver`"]
lib = CCompilerDriver`CreateLibrary["
#include <vector>
#include <algorithm>
#include <WolframCompileLibrary.h>

using namespace std;

vector<mint> sort_add2(vector<mint> vec)
{
    sort(vec.begin(), vec.end());
    for (mint &i : vec)
    {
        i = i + 2;
    }
    return vec;
}

EXTERN_C DLLEXPORT MTensor sort_add2_fc(MTensor inmma)
{
    auto pinmma = (mint*)inmma->data;//or MTensor_getIntegerDataMacro(inmma);
    vector<mint> inc(pinmma, pinmma + inmma->nelems);
    vector<mint> outc = sort_add2(inc);
    copy(outc.begin(),outc.end(), pinmma);
    (inmma->refcount)++;
    return inmma;
}", "sortAdd2File", Language -> "C++"];

dec = LibraryFunctionDeclaration[fun -> "sort_add2_fc", 
   lib, {"PackedArray"::["MachineInteger", 1]} -> 
    "PackedArray"::["MachineInteger", 1]];
fc = FunctionCompile[dec, 
  Function[Typed[n, "PackedArray"::["MachineInteger", 1]], fun[n]]]

Then we can use it now:

fc[{2, 3, 1}]

{3, 4, 5}

---

Performace with Ben Izd's fn:

Needs["GeneralUtilities`"]
GeneralUtilities`BenchmarkPlot[<|"Sort[lis]+2" -> (Sort[#] + 2 &), 
  "yode" -> fc, "Ben Izd" -> fn|>, RandomInteger[100, #] &, 
 Range[100, 1000000, 1000], "IncludeFits" -> True]

---

Add another case to calculate the mean for a ragged list here. As the direction of the author, we can greatly boost the speed by a factor of 5 than the uncompiled solution. It is a significant dig, so I share it here for the community.

The huge trouble is the c++ cannot receive the "ListVector"["PackedArray"["Real64", 1]], we can see the compile string with FunctionCompileExportString for the format:

FunctionCompileExportString[
 Typed[x, "ListVector"["PackedArray"["Real64", 1]]] |-> x]

line6~9

We have built a struct receive the parameter from MMA in c++ like:

struct ManagedArray {
    MTensor* data;
    mint length;
};

struct FixedArray{
    ManagedArray data;
    mint refcount;
};

struct RaggedTensor{
    FixedArray* data;
    mint refcount;
};

So the memory for both the parameters and the arguments are aligned. Here is the complete code:

Needs["CCompilerDriver`"]
lib=CCompilerDriver`CreateLibrary["
#include<vector>
#include<numeric>
#include<WolframCompileLibrary.h>

using namespace std;

struct ManagedArray {
    MTensor* data;
    mint length;
};

struct FixedArray {
    ManagedArray data;
    mint refcount;
};

struct RaggedTensor {
    FixedArray* data;
    mint refcount;
public:
    mint size() {
        return data->data.length;
    }

    MTensor* begin() {
        return data->data.data;
    }

    MTensor* end() {
        return this->begin() + this->size();
    }
    vector<vector<mreal>> toVecs() {
        vector<vector<mreal>> vecs;
        for (auto tensor : *this) {
            auto ptensor = (mreal*)tensor->data;
            vecs.push_back(vector<mreal>(ptensor, ptensor + tensor->nelems));
        }
        return vecs;
    }
};

mreal mean(vector<vector<mreal>> vecs) {
    vector<mreal> tmp;
    for (vector<mreal> vec : vecs) {
        tmp.insert(tmp.end(), vec.begin(), vec.end());
    }
    return accumulate(tmp.begin(), tmp.end(), 0.) / tmp.size();
}

EXTERN_C DLLEXPORT mreal ragged_mean(RaggedTensor in) {
    return mean(in.toVecs());
}","Ragged_Mean",Language->"C++"];

dec=LibraryFunctionDeclaration[fun->"ragged_mean",lib,{"ListVector"::["PackedArray"::["Real64",1]]}->"Real64"];
fc=FunctionCompile[dec,Function[{Typed[lis,"ListVector"::["PackedArray"::["Real64",1]]]},fun[lis]]]

We can use it now:

raglist = N[{{5, 9, 3}, {7, 9}, {9, 11, 8, 33}}];
fc[Flatten[raglist], Length /@ raglist]

10.4444

Performace:

lis=DeleteCases[Table[ResourceFunction["RandomSplit"][RandomReal[255, 1000], 10], 
    100000], {}, {2}];

RepeatedTiming[fc /@ lis;]
RepeatedTiming[Mean@*Flatten /@ lis;]

{0.673275, Null}

{2.80634, Null}