19
$\begingroup$

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

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

$\endgroup$
5
  • 1
    $\begingroup$ I don't know the answer, but as a side remark you can see when anything gets unpacked by first switching on the message like so: On[Developer`FromPackedArray::punpack1]. Then evec = Eigenvectors[mat]; gives a warning Developer`FromPackedArray::punpack1: Unpacking array with dimensions {300}. >> $\endgroup$
    – acl
    Nov 2, 2013 at 21:39
  • $\begingroup$ @acl Thanks, I know this feature! On["Packing"] is more convenient. $\endgroup$
    – ybeltukov
    Nov 2, 2013 at 21:47
  • $\begingroup$ Would it be an option to call the MKL by yourself? The Intel Compiler is free as academic version for Linux. If it is really numerical calculation of EigenVectors only, it could be worth it to call the LAPACK routing by yourself. $\endgroup$
    – halirutan
    Nov 2, 2013 at 22:36
  • $\begingroup$ @halirutan Thank you for the idea. Several years ago I used this library. However, it would be great to do whole process (preparation of the matrix, solving the eigensystem, and further analysis) in Mathematica. Such very big matrices are rare but even medium-size (~1/10 of RAM) matrices can make a system go into swap and this is annoying. I can't visit this site at this moments :) $\endgroup$
    – ybeltukov
    Nov 2, 2013 at 23:05
  • 1
    $\begingroup$ I agree this unpacking is not a nice thing to happen. We'll look into it. $\endgroup$ Nov 4, 2013 at 16:34

1 Answer 1

12
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Thank you for the great answer! I did not imagine that it is so easy to connect Mathematica and MKL. I hope this problem will be solved in the next versions of Mathematica. $\endgroup$
    – ybeltukov
    Nov 4, 2013 at 17:16
  • 2
    $\begingroup$ @ybeltukov I think it will be addressed. $\endgroup$ Nov 4, 2013 at 18:48
  • 1
    $\begingroup$ This seems like the best option if one has MKL. If not, perhaps LinearAlgebra`LAPACK`GEEV could be useful? The problem is that I'm not sure which 9 of the 13 arguments for *GEEV are actually required, or in what form and order. $\endgroup$ Nov 6, 2013 at 20:21
  • $\begingroup$ I just realized that actually you don't use MKL! You use only definitions, but shared library comes from Mathematica itself. You can put this definitions instead #include <mkl_lapacke.h> and everething will work. $\endgroup$
    – ybeltukov
    Sep 4, 2014 at 20:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.