Casting complex to real without data copy in MATLAB R2018a and newer

Question

Since MATLAB R2018a, complex-valued matrices are stored internally as a single data block, with the real and imaginary component of each matrix element stored next to each other -- they call this "interleaved complex". (Previously such matrices had two data blocks, one for all real components, one for all imaginary components -- "separate complex".)

I figure, since the storage now allows for it, that it should be possible to cast a complex-valued array to a real-valued array with twice as many elements, without copying the data.

MATLAB has a function typecast, which casts an array to a different type without copying data. It can be used, for example, to cast an array with 16 8-bit values to an array with 2 double floats. It does this without copying the data, the bit pattern is re-interpreted as the new type.

Sadly, this function does not work at all on complex-valued arrays.

I'm looking to replicate this code:

A = fftn(randn(40,60,20)); % some random complex-valued array
assert(~isreal(A))

sz = size(A);
B = reshape(A,1,[]);        % make into a vector
B = cat(1,real(B),imag(B)); % interleave real and imaginary values
B = reshape(B,[2,sz]);      % reshape back to original shape, with a new first dimension
assert(isreal(B))

The matrices A and B have (in R2018a and newer) the exact same data, in exactly the same order. However, to get to B we had to copy the data twice.

I tried creating a MEX-file that does this, but I don't see how to create a new array that references the data in the input array. This MEX-file works, but causes MATLAB to crash when clearing variables, because there are two arrays that reference the same data without them realizing that they share data (i.e. the reference count is not incremented).

// Build with:
//    mex -R2018a typecast_complextoreal.cpp

#include <mex.h>

#if MX_HAS_INTERLEAVED_COMPLEX==0
#error "This MEX-file must be compiled with the -R2018a flag"
#endif

#include <vector>

void mexFunction(int /*nlhs*/, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {

   // Validate input
   if(nrhs != 1) {
      mexErrMsgTxt("One input argument expected");
   }
   if(!mxIsDouble(prhs[0]) && !mxIsSingle(prhs[0])) {
    mexErrMsgTxt("Only floating-point arrays are supported");
   }

   // Get input array sizes
   mwSize nDims = mxGetNumberOfDimensions(prhs[0]);
   mwSize const* inSizes = mxGetDimensions(prhs[0]);

   // Create a 0x0 output matrix of the same type, but real-valued
   std::vector<mwSize> outSizes(nDims + 1, 0);
   plhs[0] = mxCreateNumericMatrix(0, 0, mxGetClassID(prhs[0]), mxREAL);

   // Set the output array data pointer to the input array's
   // NOTE! This is illegal, and causes MATLAB to crash when freeing both
   // input and output arrays, because it tries to free the same data
   // twice
   mxSetData(plhs[0], mxGetData(prhs[0]));

   // Set the output array sizes
   outSizes[0] = mxIsComplex(prhs[0]) ? 2 : 1;
   for(size_t ii = 0; ii < nDims; ++ii) {
      outSizes[ii + 1] = inSizes[ii];
   }
   mxSetDimensions(plhs[0], outSizes.data(), outSizes.size());
}

I'd love to hear of any ideas on how to proceed from here. I don't necessarily need to fix the MEX-file, if the solution is purely MATLAB code, so much the better.

Answer 1

See this FEX submission, which can do the complex --> 2 reals reinterpretation without copying the data (it can even point to interior contiguous sub-sections of the data without a copy):

https://www.mathworks.com/matlabcentral/fileexchange/65842-sharedchild-creates-a-shared-data-copy-of-a-contiguous-subsection-of-an-existing-variable

If you are simply reading and writing interleaved complex data files in R2018a and later, see this FEX submission:

https://www.mathworks.com/matlabcentral/fileexchange/77530-freadcomplex-and-fwritecomplex

Answer 2

This question reminded me of some blog posts regarding in-place memory editing through MEX, about which the following was said:

[M]odifying the original data directly is both discouraged and not officially supported. Doing it incorrectly can easily crash Matlab. ^ref

This is, in the best case, an unmaintainable mess. ^ref

Having said that, I don't have a solution for you, but I might be able to suggest a workaround.

Seeing how MATLAB allows us to call python libraries, we can perform these sort of manipulations within python, and only bring the data back into MATLAB when we reach a stage where proceeding in python is impossible or undesired. Depending on when this stage might be, this idea could be either a valid approach, or a completely useless suggestion.

Take a look at the example below, it should be pretty self-explanatory:

np = py.importlib.import_module('numpy');
sp = py.importlib.import_module('scipy.fftpack');

% Create a double array in python:
arrC = sp.fftn(np.random.rand(uint8(4), uint8(3), uint8(2)));
%{
arrC = 
  Python ndarray with properties:

           T: [1×1 py.numpy.ndarray]
        base: [1×1 py.NoneType]
      ctypes: [1×1 py.numpy.core._internal._ctypes]
        data: [1×4 py.memoryview]
       dtype: [1×1 py.numpy.dtype]
       flags: [1×1 py.numpy.flagsobj]
        flat: [1×1 py.numpy.flatiter]
        imag: [1×1 py.numpy.ndarray]
    itemsize: [1×1 py.int]
      nbytes: [1×1 py.int]
        ndim: [1×1 py.int]
        real: [1×1 py.numpy.ndarray]
       shape: [1×3 py.tuple]
        size: [1×1 py.int]
     strides: [1×3 py.tuple]
    [[[ 13.99586491+0.j           0.70305071+0.j        ]
      [ -1.33719563-1.3820106j   -0.74083670+0.25893033j]
      [ -1.33719563+1.3820106j   -0.74083670-0.25893033j]]

     [[ -0.43914391+0.8336674j    0.08835445-0.50821244j]
      [  1.07089829-0.35245746j   0.44890850-0.9650458j ]
      [  2.09813180+1.34942678j  -1.20877832+0.71191772j]]

     [[ -2.93525342+0.j          -0.69644042+0.j        ]
      [  0.16165913-1.29739125j  -0.84443177+0.26884365j]
      [  0.16165913+1.29739125j  -0.84443177-0.26884365j]]

     [[ -0.43914391-0.8336674j    0.08835445+0.50821244j]
      [  2.09813180-1.34942678j  -1.20877832-0.71191772j]
      [  1.07089829+0.35245746j   0.44890850+0.9650458j ]]]
%}

% Make sure that python sees it as a "complex double" (aka complex128)
assert( isequal(arrC.dtype, np.dtype(np.complex128)) );

% Return a (real) double view:
arrR = arrC.view(np.float64);
%{
arrR = 
  Python ndarray with properties:

           T: [1×1 py.numpy.ndarray]
        base: [1×1 py.numpy.ndarray]
      ctypes: [1×1 py.numpy.core._internal._ctypes]
        data: [1×4 py.memoryview]
       dtype: [1×1 py.numpy.dtype]
       flags: [1×1 py.numpy.flagsobj]
        flat: [1×1 py.numpy.flatiter]
        imag: [1×1 py.numpy.ndarray]
    itemsize: [1×1 py.int]
      nbytes: [1×1 py.int]
        ndim: [1×1 py.int]
        real: [1×1 py.numpy.ndarray]
       shape: [1×3 py.tuple]
        size: [1×1 py.int]
     strides: [1×3 py.tuple]
    [[[ 13.99586491   0.           0.70305071   0.        ]
      [ -1.33719563  -1.3820106   -0.7408367    0.25893033]
      [ -1.33719563   1.3820106   -0.7408367   -0.25893033]]

     [[ -0.43914391   0.8336674    0.08835445  -0.50821244]
      [  1.07089829  -0.35245746   0.4489085   -0.9650458 ]
      [  2.0981318    1.34942678  -1.20877832   0.71191772]]

     [[ -2.93525342   0.          -0.69644042   0.        ]
      [  0.16165913  -1.29739125  -0.84443177   0.26884365]
      [  0.16165913   1.29739125  -0.84443177  -0.26884365]]

     [[ -0.43914391  -0.8336674    0.08835445   0.50821244]
      [  2.0981318   -1.34942678  -1.20877832  -0.71191772]
      [  1.07089829   0.35245746   0.4489085    0.9650458 ]]]
%}

% Do something else with it in python
...

% Bring data to MATLAB:
sz = cellfun(@int64, cell(arrR.shape));
B = permute(reshape(double(py.array.array('d', arrR.flatten('C').tolist())),sz),[2,1,3]);

The last stage might be inefficient enough to negate the performance/memory gains by not copying the data previously, so it would probably be a good idea to keep it for the very end, and reduce the data as much as possible beforehand.

---

While experimenting with this, I came to realize that there's much difficulty¹ involved in transferring non-scalar complex data from MATLAB to Python and back (just compare the MATLAB output for np.asarray(1+1i) vs np.asarray([1+1i, 1+1i])). Possibly, the reason for this is the limited support for complex non-ndarray arrays in python.

If you (as opposed to me) know what to do with memoryview objects (i.e. the content of the data field of the ndarray objects) - you could get a pointer and perhaps pass this on to C to get some useful results.

---

¹ It's possible, but you have to num2cell your data so that it can get transferred as a python list (and vice versa). The situation can also be improved by editing some MATLAB files.