Hello,

I'm pretty new to eigen and am not a natural-born programmer, either (I am a physicist). However, I am doing a lot of image reconstruction in matlab and MRI sequence development in C++. Sometimes it is very useful to share code between these languages - that's why I'm working with matlab's mex a lot (using the gcc compiler).

When I found out about eigen, I was pretty impressed with its feature set and ease of use (it' really not that much harder than matlab). My problem right now is the interface between matlab and eigen/C++. It is pretty easy to copy data from a matlab-array to an "eigen"-matrix (no pun intended...) using matlab's mex.h library. Since I'm working with fairly big data sets, mapping the data to the MatrixXd object is of course much better memory wise. I already found out how to do that with real data:

Code: Select all

double* pData = mxGetPr(prhs[0]);
int m   = mxGetM(prhs[0]);
int n   = mxGetN(prhs[0]);
MatrixXd data = Map<MatrixXd>(pData, m, n);


(using some of matlab's mex functions, provided by mex.h). This works exactly as expected and I'm quite happy with that.

So, here comes the question: Is it at all possible to map arrays to a complex MatrixXcd object? Matlab basically provides two pointers for each complex array, one for the real part, one for the imaginary part. I need to map these two now to MatrixXcd.

Any help is appreciated!

Somebody please correct me if I'm wrong, but here's my understanding:

Eigen complex matrices are Matrix<complex>, but it looks like the MATLAB storage format is complex<array>. That is, for Eigen
matrices, the real and complex elements are interleaved in memory, whereas for MATLAB, there is a contiguous chunk for the reals followed by a contiguous chunk for the imaginaries.

You could read from the complex MATLAB arrays:

Code: Select all

double *preal=mxGetPr(prhs[0]), *pimag=mxGetPi(prhs[0]);

int m   = mxGetM(prhs[0]), n   = mxGetN(prhs[0]);

MatrixXd realPart = Map<MatrixXd>(preal, m, n),
                 imagPart = Map<MatrixXd>(pimag, m, n);

MatrixXcd data = realPart + complex<double>(0,1.) * imagPart;


I'm not sure if data would be an interleaved copy of the MATLAB array ... if it was, you could do some math and then write back to the original arrays through

Code: Select all

realPart = data.real();
imagPart = data.imag();


note that when you do:

MatrixXd mat = Map<MatrixXd>(pData, m, n);

the MatLab data are copied into mat. You can avoid this copy by naming the Map object:

Map<MatrixXd> mat(pData, m, n);

You can use mat just like a MatrixXd with the difference that all your own function taking Eigen's objects as argument are correctly templated, i.e., do not declare:

void foo(const MatrixXd& a);

but:

template<typename A>
void foo(const MatrixBase<A>& a);

Now regarding complexes, well, if the imaginary and real parts are stored in different array and not interleaved then there is no alternative than copying the data to a MatrixXcd as follow:

MatrixXcd mat(m,n);
mat.real() = Map<MatrixXd>(realData,m,n);
mat.imag() = Map<MatrixXd>(imagData,m,n);

I'm surprised that MatLab uses this kind of storage because it not BLAS/LAPACK compatible, and as far as I know MatLab relies on BLAS/LAPACK, so I'm puzzled...

Thanks a lot for your quick replies! My code runs already noticeably faster

Unfortunately, I don't know much about matlab's internals - the mex.h file is also not very helpful. Here is what I know:
The mex-interface provides a struct called mxArray (it includes the real part, imaginary part, and some information like type, array size...). With the functions mxGetPr and mxGetPi I get access to the real and imaginary part, respectively (and can use them like any other C-style array). So i would assume that they are stored as separate entities?

indeed, it seems that's right:
http://www.mathworks.cn/access/helpdesk ... tml#f12246

so you can either copy the buffers to a MatrixXcd as I showed you, or map the real and imaginary parts and do the complex arithmetic by hand, e.g., a matrix-matrix complex product requires 4 real-products... it depends on what you have to do.