Unpacked eigenvectors of complex matrices

Question

Edit: This bug is fixed in Mathematica 10.0.0

I want to calculate eigenvectors of a Hermitian matrix. For example

$HistoryLength = 0;
n = 3000;
mat = # + ConjugateTranspose[#] &@RandomComplex[{0, 1 + I}, {n, n}];
evec = Eigenvectors[mat];

However, evec is unpacked!

evec // Developer`PackedArrayQ

False

System monitor shows

It seems that the result was unpacked in the end of evaluation. I think the problem is that the last element of each eigenvector is real instead of complex

evec[[All, -1]]

-0.100012, -0.059244, -0.0723702, ...

evec[[All, -2]]

-0.0994738 - 0.00136552 I, -0.0527357 - 0.0532419 I, 0.158174 + 0.0341157 I, ...

With real matrices everything is fine

n = 5000;
mat = # + Transpose[#] &@RandomReal[1, {n, n}];
evec = Eigenvectors[mat];
evec // Developer`PackedArrayQ

True

My current workaround consist in packing the result immediately after Eigenvectors

evec = Developer`ToPackedArray[0.0 I + Eigenvectors[mat]];

Now evec is packed but there is a very big peak in the memory usage

It is not a problem for a small matrices, but if matrix mat is ~30GB then I need ~300GB memory. It is quite a lot even for the swap space!

It seems that this is a very simple question, but I could not find anything.

P.S. I use Mathematica 9.0.1 and Linux. Mathematica 8 has the same problem.

Update

halirutan's answer solves this problem completely! Now memory profile is much better

Some notes:

You need 11 version of MKL. In Gentoo it is available in science overlay. I compile the library with the command

gcc -shared -fPIC -o PackedEigenvectors.so -I /usr/local/Wolfram/Mathematica/9.0/SystemFiles/IncludeFiles/C/ -I /opt/intel/composerxe-2013.0.080/mkl/include/ PackedEigenvectors.c

Helpful links: 1, 2

Answer 1

Let me put my comment into an answer, because I think we might have misunderstood each other. You answered in the comment

However, it would be great to do whole process (preparation of the matrix, solving the eigensystem, and further analysis) in Mathematica.

That exactly was my idea. You only write some lines of C-Code which are compiled into a WolframLibrary function. This can directly be called from within Mathematica. So the only thing you do is, you substitute the call of Eigenvectors. For the following, please note that I haven't used the MKL in a while and I basically just copied an example which looked promising.

Here is the complete Library wrapper code for the MKL function which calculates the eigensystem of a complex hermitian matrix. The core of the wrapper is the following function

DLLEXPORT int zheev(
        WolframLibraryData libData, 
        mint argc, 
        MArgument *args,
        MArgument res) {
    MTensor input = MArgument_getMTensor(args[0]);
    MKL_Complex16 *data = (MKL_Complex16*) libData->MTensor_getComplexData(input);
    const mint *dims = libData->MTensor_getDimensions(input);

    MKL_INT n = dims[0], info, lda = n;

    double *w = (double*) malloc(n*sizeof(double));

    info = LAPACKE_zheev(LAPACK_ROW_MAJOR, 'V', 'L', n, data, lda, w);

    if (info > 0) {
        free(w);
        return LIBRARY_NUMERICAL_ERROR;
    }

    MArgument_setMTensor(res, input);
    return LIBRARY_NO_ERROR;

}

This library function expects a complex matrix as input whose numeric values are stored in data. Everything depends here on the fact that the WolframLibrary stores its mcomplex in the same way as the MKL_Complex16. In my case here this means to check whether the memory layout for the following types are identical

typedef struct {mreal ri[2];} mcomplex;

typedef struct _MKL_Complex16 {
    double real;
    double imag;
} MKL_Complex16;

mreal is typedef'd as double here. If this is the case then chances are very good that you can

• extract the tensor data from the input matrix (getting a C array of complexes)
• feed this array directly to the MKL routine which stores the eigenvectors by overwriting exactly this array
• return your input tensor back to Mathematica

If the C part works, then installing and calling the WolframLibrary is as easy as the next few lines of code

AppendTo[$LibraryPath, "/path/where/your/lib/is"];
fun = LibraryFunctionLoad["libCallMKL", "zheev", {{_Complex, 2}}, {_Complex, 2}];

On["Packing"]
n = 3000;
mat = # + ConjugateTranspose[#] &@RandomComplex[1, {n, n}];
evec = fun[mat];

In my tests, this took fewer memory as the Eigenvectors call and it did not unpack the array in Mathematica.