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
On[Developer`FromPackedArray::punpack1]
. Thenevec = Eigenvectors[mat];
gives a warningDeveloper`FromPackedArray::punpack1: Unpacking array with dimensions {300}. >>
$\endgroup$On["Packing"]
is more convenient. $\endgroup$