contiguous memory for complex arrays in mex

Is it possible to allocate a complex array in a mex file and have the real and imaginary parts be contiguous with each other? I'm working with very large matrices and need to conserve memory and would like to not have to make a copy of the input matrices to get them into FORTRAN format. If the memory were contiguous, I could just rearrange the data into interleaved real and imaginary. I've already been able to pass in a real valued array from matlab and have LAPACK treat it as a complex (it just assumes the real and imaginary are interleaved), but this isn't exactly an elegant solution since I still have to write the data to a file and read it back into complex form. At that point it becomes useless to use MATLAB at all and I might as well just write a standalone program in C or FORTRAN. Any ideas if this can be done or am I out of luck?

Answer 1

[James Tursa]

Yes and No. It is easily possible to physically do this, but the clearing of such a variable will likely crash MATLAB if you are not very careful. The outline is as follows.

In the mex routine create a double array that is sized to hold all of your intended complex data. E.g., if you want a 1x5 complex array then create a double array that is 1x10. Then manually use mxSetN to change the column size to 5. Then pick the midpoint of the real data area and use mxSetPi to set the imaginary data pointer to that address. That will create a complex mxArray variable with the real and imaginary parts contiguous.

This will work fine up until the point that MATLAB tries to clear this variable from memory. Then it will try to free the address you just set with the mxSetPi function. Since this is not a valid address to free, the results are unpredictable and MATLAB will likely crash.

The only robust solution I can think of is to manually create a shared data copy of this variable, make it persistent, and keep this shared data copy inside the mex routine so MATLAB cannot try to clear it. Then use a mexAtExit function to manually reverse the process, use mxSetPi with a NULL, and then clear the variable. But I will tell you right now that this is very tricky to accomplish and requires that you dig under the hood of the mxArray structure to do some unofficial manipulation. You need to be able to detect when the mxArray is shared with other variables inside your mex routine, etc etc. Not recommended for the faint of heart and not easy to accomplish even for a seasoned programmer IMO.

You might be able to get away with it if you can keep all the operations inside the mex routine so that you can carefully manage all of the operations on the variable to ensure there is no shared data copy of it floating around in memory elsewhere, and then manually fix & clear it yourself inside the mex routine.

What operations will you be doing with this variable?

[Steve]

Thanks James, that is exactly what I thought/was afraid of. I think it's probably not worth the effort to do this in MATLAB and I should just make a stand alone program that can write out a file with the result I need. I plan on solving a linear system of equations, but using the packed-symmetric form of the matrix. This lets me get about 40% more unknowns into memory then making the full explicit matrix in MATLAB.

Just to satisfy my curiosity, do you have any idea why the following code doesn't work? It instantly crashes MATLAB, not after I clear the variable. Oddly enough it works with ELEMENTS=4, but not any bigger.

#define ELEMENTS 6

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

double *dynamicData; // pointer to dynamic data //

mwSize index;

const double data[] = {2.1, 3.4, 2.3, 2.45, 1.1, 2.2};

// Create a local array and load data //

dynamicData = mxMalloc(ELEMENTS * sizeof(double));

for ( index = 0; index < ELEMENTS; index++ ) {

dynamicData[index] = data[index];

}

// Create a 0-by-0 mxArray //

plhs[0] = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxREAL);

// Put the C array into the mxArray and define its dimensions //

mxSetPr(plhs[0], dynamicData);

mxSetM(plhs[0], 1);

mxSetN(plhs[0], ELEMENTS/2);

mxSetPi(plhs[0], &dynamicData[ELEMENTS/2]);

return;

}

[James Tursa]

This crashes MATLAB because the mxSetPr, mxSetPi, mxSetEtc... functions all check to see that the address passed to them originally came from the MATLAB memory manager (i.e., the result of a mxMalloc, mxCalloc, etc call). If not, then it crashes intentionally (with an assertion fault or similar if I remember correctly). Earlier versions of MATLAB did not make this check, but later versions of MATLAB do. The address &dynamicData[0] is known to the API so it would work in such calls, but &dynamicData[ELEMENTS/2] is not. That is why you have to dig under the hood of the mxArray to manually set the value of pi by bypassing the checks in the API functions. But then you have the clearing problem downstream (which is what these API checks are trying to avoid). It is not a bad thing for the API to do since it forces a crash at the point the problem is created rather than downstream where it would be hard for the programmer to understand where the root fault is.

[James Tursa]

P.S. I am surprised that it works for ELEMENTS=4. I don't understand that one.

[James Tursa]

P.P.S. This problem intrigued me so I went ahead and coded it up. I will write up some documentation and post it to the FEX in the next day or two. My routine will basically turn any numeric real variable (with an even number of elements) into a "contiguous" complex variable in-place, or vice-versa. Keeps a shared data copy of it off to the side so that MATLAB will not attempt to clear it. The routine will clean up everything once the user is done with the variables.

[Sean Little]

James, did you ever post that code? It sounds like exactly what I need, but I don't see it on the file exchange.

Answer 2

[the cyclist]

This page has info that looks like it might be helpful for what you are trying to do:

http://www.mathworks.com/help/techdoc/matlab_external/br_2m24-1.html#f32627

In particular, it mentions the fort2mat() and mat2fort() functions.

[Steve]

Thanks cyclist. I've used those functions before, but they have the problem of having to create a copy of the data before they can be used. If I wanted to do an inversion, I'd have to create a copy of the input array using mat2fort, use lapack to invert it, and create another copy using fort2mat to return the answer. I'd need 3x the memory I would if I could just force the pointer to the imaginary part of the matrix to come immediately after the real part and do the inversion in place.