3
$\begingroup$

I want to use C++ Eigen library in Mathematica Librarylink function.

Currently, I have no problem dealing with real arrays.

But I don't know what is the best way to deal with complex array in Librarylink in terms of C++.

So it is a C++ question in essence, but it is also a librarylink question.

In mma librarylink, the header file WolframLibrary.h provide type mcomplex, it is actually a structure which has length 2 real array representing real and imaginary part respectively, it is defined as

typedef struct {mreal ri[2];} mcomplex;

I don't know whether it is standard or not, but it seems awkward, and more importantly it is not directly compatible with Eigen.

While in C++, there is standard header <complex>. After including this header, we can define complex number simply as

std::complex<double> cc(1.0,2.0);

Eigen library directly compatible with this complex type. For example, we could do the following

#include <iostream>
#include <complex>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;
int main(){
    complex<double> *cc;    //define complex double pointer
    Matrix2cd m=Matrix2cd::Random();   //Eigen 2x2 complex matrix object
    cc=m.data();    //m.data() return a C array pointer and assign it to cc
    for(int i=0;i<4;i++)
    cout<<cc[i]<<endl;   //output cc array, the data that cc points to
                         //is the same as Eigen object m
}

Now I give an example of using C++ Eigen in librarylink function, served for testing and reference.

libnewsource = "#include<Eigen/Dense>
  #include<Eigen/LU>
  #include \"WolframLibrary.h\"
  DLLEXPORT mint WolframLibrary_getVersion(){return \
WolframLibraryVersion;}
  DLLEXPORT int WolframLibrary_initialize( WolframLibraryData \
libData) {return 0;}
  DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData \
libData) {}

  EXTERN_C DLLEXPORT int eigeninverse(WolframLibraryData libData, \
mint Argc, MArgument *Args, MArgument Res) {

  using namespace Eigen;

  int err; // error code
  //input-------------
  MTensor m1;
  double *ptm1;
  const mint *dim;
  //output------------
  MTensor outT;
  double *out;

  //input------------------
  m1 = MArgument_getMTensor(Args[0]);
  ptm1=libData->MTensor_getRealData(m1);
  dim=libData->MTensor_getDimensions(m1);   
  //output--------------------
  err = libData->MTensor_new(MType_Real, 2, dim, &outT);
  out=libData->MTensor_getRealData(outT);

  Map<MatrixXd> eigenm(ptm1,dim[0],dim[1]);  //use Map to convert \
array pointed by ptm1 to Eigen matrix type

  Map<MatrixXd>(out,dim[0],dim[1])=eigenm.inverse();  \
//use Map to convert Eigen matrix into out array

  MArgument_setMTensor(Res,outT);
  return LIBRARY_NO_ERROR;
  }";

Now we can load it with below code

LibraryFunctionUnload[eigeninverse];
libnew = CreateLibrary[libnewsource, "eigeninverse","CompileOptions"->"-x c++"];
eigeninverse = 
  LibraryFunctionLoad[libnew, "eigeninverse", {{Real, 2}}, {Real, 2}];

About compiler setting:

For mingw user: make sure to set compiler as x86_ 64-w64-mingw32-g++.exe, not x86_64-w64-mingw32-gcc.exe. This is what I use in init.m file

CCompilerDriver`$CCompiler=
  {"Compiler"->GenericCCompiler,
   "CompilerInstallation"->"C:/mingw-w64/mingw64",
   "CompilerName"->"x86_64-w64-mingw32-g++.exe",
   "CompileOptions"->"-O2 -x c"}

To notice there is "CompileOptions"->"-O2 -x c". -x c is g++ option to treat file as c. I found without this option, normal Compile will give errors like

"The function compiledFunction0 was not loaded from the file "somepath/.../compiledFunction0.dll "

And to create librarylink function that written in C++, you should add CompileOptions"->"-x c++ to CreateLibrary to threat file as C++.


My Question

What is the best way to transform the above code to support complex matrix. The problem mainly focus on how to convert efficiently between Mathematica mcomplex and C++ standard complex<double>. And I personally feel that dealing with real and imaginary part separately kind of awkward. I don't know whether it is possible or not that we modify WolframLibrary.h to make it support C++ standard?

$\endgroup$
  • $\begingroup$ If you're going to do a lot of LibraryLink stuff, you should really start putting the code in separate files. Don't use strings! I believe since C++11 the layout of std::complex<double> is guaranteed. It is identical to mcomplex. Thus you can reinterpret_cast one to the other when needed. $\endgroup$ – Szabolcs Jan 12 '16 at 11:33
  • $\begingroup$ @Szabolcs Hi, szabolcs. OK, I will if the code getting larger. But what do you mean by "identical to mcomplex", and how to use reinterpret_cast $\endgroup$ – matheorem Jan 12 '16 at 11:38
  • $\begingroup$ @Szabolcs I don't understand. In C99, the complex number is defined by macro gnu.org/software/libc/manual/html_node/Complex-Numbers.html It is not the same as mathematica's mcomplex, is mcomplex also a standard? $\endgroup$ – matheorem Jan 12 '16 at 13:26
  • 1
    $\begingroup$ @matheorem reinterpret_cast is the C++ operator that is, in many ways, most similar to old-style C casting. The code reinterpret_cast<new_type>(value) casts value into the type new_type much as (new_type) value would have in C. By "identical to mcomplex" Szabolcs just means that there isn't a difference in how the two are stored in memory, only in how they are handled. Both are stored of a pair of double-values in memory; because of this, just reinterpreting one as the other will work (it wouldn't work if, for example, one of them used float and the other used double). $\endgroup$ – nben Jan 12 '16 at 22:03
  • $\begingroup$ @user21382 Would you like to make an answer? $\endgroup$ – matheorem Jan 13 '16 at 2:44
5
$\begingroup$

LibraryLink's complex type, mcomplex, is defined as two contiguous double values.

In C++, std::complex<double> has exactly the same layout. This is guaranteed since C++11, but should hold in most other cases too in practice.

This means that if you get a complex array from Mathematica,

mcomplex *arr = MTensor_getComplexData(t);

then you can reinterpret it as an array of std::complex<double> values because it has the exact same layout in memory as those:

std::complex<double> *ca = reinterpret_cast< std::complex<double> * >(arr);

I won't comment on Eigen because I am not familiar with it.

$\endgroup$
  • $\begingroup$ Thank you so much, Szabolcs. You've helped me a lot in my learning of librarylink! In terms of reintrepret_cast, I suppose it should not have overhead, just reinteprete of data, right? $\endgroup$ – matheorem Jan 14 '16 at 7:02
  • $\begingroup$ also I have pasted a complex version of my original real version of Inverse. I find it the data exchange part is really cumbersome. I am not sure whether I have written something unnecessary or not, could you please have a look at it? $\endgroup$ – matheorem Jan 14 '16 at 7:05
  • 1
    $\begingroup$ And what happens if Mathematica changes the mcomplex type ;-) $\endgroup$ – user21 Jan 14 '16 at 7:23
  • $\begingroup$ @user21 Will it? I can't find any documentation on how to access the real and imaginary parts of an mcomplex at all, but without being able to access them it isn't really possible to work with complexes ... so it's not unreasonable to resort to looking at the header and assume that users are meant to access the struct directly. $\endgroup$ – Szabolcs Jan 14 '16 at 8:38
  • $\begingroup$ @user21 Hm, now I found mcreal and mcimag in demo.c and in WolframLibrary.h. $\endgroup$ – Szabolcs Jan 14 '16 at 8:39
1
$\begingroup$

Thanks to Szabolcs for pointing out the right way.

So I rewrite the real Inverse version into complex version, and pasted here for reference.

libnewsource="#include<complex>
#include<Eigen/Dense>
#include<Eigen/LU>
#include \"WolframLibrary.h\"
DLLEXPORT mint WolframLibrary_getVersion(){return WolframLibraryVersion;}
DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {return 0;}
DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {}

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

using namespace Eigen;
using namespace std;

int err; // error code

//input data-------------
MTensor m1=MArgument_getMTensor(Args[0]);
mcomplex *arr=libData->MTensor_getComplexData(m1);
complex<double> *ca = reinterpret_cast<complex<double>*>(arr);
const mint *dim=libData->MTensor_getDimensions(m1);

//output data--------------------
MTensor outT;
err = libData->MTensor_new(MType_Complex, 2, dim, &outT);
mcomplex *out;
out=libData->MTensor_getComplexData(outT);
complex<double> *outcom=reinterpret_cast<complex<double>*>(out);

//-----------------code part-------------------------

Map<MatrixXcd> eigenm(ca,dim[0],dim[1]);
Map<MatrixXcd>(outcom,dim[0],dim[1])=eigenm.inverse(); 

//-----------------code part-------------------------

MArgument_setMTensor(Res,outT);
return LIBRARY_NO_ERROR;
}";

and load it as below, don't foget to put -std=c++11

LibraryFunctionUnload[eigeninverse];
libnew = CreateLibrary[libnewsource, "eigeninverse", 
   "CompileOptions" -> "-x c++ -std=c++11"];
eigeninverse = 
  LibraryFunctionLoad[libnew, 
   "eigeninverse", {{Complex, 2}}, {Complex, 2}];

You can test it with

A = RandomComplex[{0. + 0. I, 1. + 1. I}, {3, 3}];
eigeninverse[A] - Inverse[A] // Chop
(*{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}*)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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