What is the best way to work with complex number array in librarylink function in terms of C++ and Eigen?

Question

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?

Answer 1

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.

Answer 2

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}}*)