How to efficiently extract real/imaginary parts of complex matrix in Eigen3 library?

Question

I have some complex, dense vectors/matrices in the Eigen3 library, and I want to extract the real and imaginary parts into separate arrays. In Matlab, I could do something like

cplxFoo = [1, 1i; -1i -1]
re = real(cplxFoo)
im = imag(cplxFoo)

which expectantly yields

cplxFoo =
   1.0000 + 0.0000i   0.0000 + 1.0000i
   0.0000 - 1.0000i  -1.0000 + 0.0000i
re =
     1     0
     0    -1
im =
     0     1
    -1     0

Is there anything like the real() and imag() Matlab functions in Eigen3?

Right now, the only thing I know will work is something akin to

MatrixXcd cplxFoo = ...;
MatrixXd re(cplxFoo.rows(), cplxFoo.cols());
MatrixXd im(cplxFoo.rows(), cplxFoo.cols());

for(size_t j=0; j<cplxFoo.cols(); ++j) {
    for(size_t i=0; i<cplxFoo.rows(); ++i) {
        re(i, j) = cplxFoo(i,j).real();
        im(i, j) = cplxFoo(i,j).imag();
    }
}

It works, and I can put it in a function even, but then I'm stuck having to do my own loop vectorizing, unrolling, etc., and I have to make an extra copy.

What I would like to be able to do is wrap a couple of Map<MatrixXd> with appropriate strides around cplxFoo to get the real and imaginary parts. But the problem is that the elements of MatrixXcd are std::complex<double>, and I'm not sure what the layout of that is. My guess is that std::complex<T> is essentially laid out like struct {T real; T imag;}; so that real and imaginary parts are tightly packed and interleaved when you make an array of std::complex<T> (and that also seems to be the consensus at this SO question), but is that guaranteed by the C++ standard? AFAICT, a compliant C++ compiler could lay it out like struct {T imag; T real;}; (note the changed order), or something more exotic like

class {
    T radius;
    T angle;

public:
    T real() const { return radius * cos(angle); }
    T imag() const { return radius * sin(angle); }
    /* ... */
};

So, is it ok to wrap a couple of Map<MatrixXd> with appropriate strides around cplxFoo? If so, how do I set up the strides correctly?

Alternatively, is there any way to get Eigen's complex data types to use separate chunks of memory for the real and imaginary parts?

For what it's worth, the reason I need to do this is because I need to interface the Eigen library with MATLAB which can only handle separate arrays for real and imaginary parts, not interleaved in any way.

Answer 1

That's easy, just use the .real() and .imag() views:

MatrixXcd M;
MatrixXd r, i;
r = M.real();
i = M.imag();

Note that you can use M.real() into an expression without copying it into a MatrixXd.