Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
1  
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[DeveloperFromPackedArray::punpack1]. Then evec = Eigenvectors[mat];` gives a warning DeveloperFromPackedArray::punpack1: Unpacking array with dimensions {300}. >>` –  acl Nov 2 '13 at 21:39
    
@acl Thanks, I know this feature! On["Packing"] is more convenient. –  ybeltukov Nov 2 '13 at 21:47
    
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. –  halirutan Nov 2 '13 at 22:36
    
@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 :) –  ybeltukov Nov 2 '13 at 23:05
    
I agree this unpacking is not a nice thing to happen. We'll look into it. –  Daniel Lichtblau Nov 4 '13 at 16:34

1 Answer 1

up vote 11 down vote accepted

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.

share|improve this answer
    
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. –  ybeltukov Nov 4 '13 at 17:16
1  
@ybeltukov I think it will be addressed. –  Daniel Lichtblau Nov 4 '13 at 18:48
1  
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. –  Oleksandr R. Nov 6 '13 at 20:21
    
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. –  ybeltukov Sep 4 at 20:53

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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