Intel MKL LAPACKE_dsyevd with n > 32766 --> Not enough memory to allocate work array in LAPACKE_dsyevd

Question

I want to compute all eigenvalues and all eigenvectors of a real symmetric matrix using LAPACKE_dsyevd from Intel MKL (2019 Update 2).

I'm having the following method in C#:

public static class MKL
{
    public static double[,] SymmetricEig(double[,] a, out double[] w)
    {
        int n1 = a.GetLength(0);
        int n2 = a.GetLength(1);
        if (n1 != n2) throw new System.Exception("Matrix must be square");
        double[,] b = Copy(a);
        int matrix_layout = 101; // row-major arrays
        char jobz = 'V';
        char uplo = 'U';
        int n = n2;
        int lda = n;
        w = new double[n];
        _mkl.LAPACKE_dsyevd(matrix_layout, jobz, uplo, n, b, lda, w);
        return b;
    }
}

with

class _mkl
{
    [DllImport(DLLName, CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)]
    internal static extern lapack_int LAPACKE_dsyevd(
        int matrix_layout, char jobz, char uplo, lapack_int n, [In, Out] double[,] a, lapack_int lda, [In, Out] double[] w);
}

and the following test code:

    int n = 32766; // 32767 or greater --> Not enough memory to allocate work array in LAPACKE_dsyevd
    double[,] A = CreateRandomSymmetricMatrix(n);
    double[] w = new double[n];
    double[,] B = MKL.SymmetricEig(A, out w);

with

    static double[,] CreateRandomSymmetricMatrix(int n1)
    {
        double[,] m = new double[n1, n1];
        for (int i1 = 0; i1 < n1; i1++)
        {
            for (int i2 = 0; i2 <= i1; i2++)
            {
                m[i1, i2] = r.NextDouble();
                m[i2, i1] = m[i1, i2];
            }
        }
        return m;
    }

If n is greater than 32766 it fails with the following error message:

Not enough memory to allocate work array in LAPACKE_dsyevd

My PC has 128 GB of RAM which should be enough. I'm aware of <gcAllowVeryLargeObjects enabled="true" /> and I've set it to true. I'm as well aware of the 65535^2 limit of a multidimensional array in C#, see 2d-Array with more than 65535^2 elements --> Array dimensions exceeded supported range.

By the way I can compute the eigenvalue decomposition using MATLAB for matrices with n = 40000 or greater. And I think MATLAB is using as well Intel MKLin the background to compute it.

So how can I compute the eigenvalue decomposition of very large matrices (n > 40000) in C# using Intel MKL?

score 0 · Answer 1 · answered Apr 05 '19 at 09:00

0

I don't think that this is your issue. The below definition suggests that w = new double[n]; is way to small in you case.

*  LWORK   (input) INTEGER
*          The dimension of the array WORK.
*          If N <= 1,               LWORK must be at least 1.
*          If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
*          If JOBZ = 'V' and N > 1, LWORK must be at least
*                                                1 + 6*N + 2*N**2.

you really should always do a workspace query. I know, that the documentation leaves this open to the user, but it's so convenient and helps to avoid your situation. So you know, how to do the workspace query? If not hit me back real quick.

answered Apr 05 '19 at 09:00

Kaveh Vahedipour

3,412
1
14
22

So just as a quick test: What happens, if you allow w to be `w = new double[2*n]`? – Kaveh Vahedipour Apr 05 '19 at 09:13
I tried it with `n=40000` and `w = new double[2*n]` which give me still the same error. – Wollmich Apr 05 '19 at 09:26
1

well, sir this seems to be a bug, then. i just dove in the code of the netlib implementation and the workspace query is done in the proper way. just cannot tell you, if the bug is on the fortran or the lapacke_ end. :( the workspace query whould return 12GB. – Kaveh Vahedipour Apr 05 '19 at 09:28
By the way: should I try to use as well `syev` or `syevr`? – Wollmich Apr 05 '19 at 09:30
What do you need? Eigenvalues and eigenvectors? What hardware (CPU)? – Kaveh Vahedipour Apr 05 '19 at 09:33
I need all Eigenvalues and Eigenvectors. My input matrix is symmetric and positive semidefinite. My PC is a HP Z440 (Intel Xeon E5-1680v4 with128GB RAM) – Wollmich Apr 05 '19 at 09:37
1

http://www.cs.utexas.edu/users/inderjit/public_papers/invit_fail.pdf describes why `dsyevr` should be the choice. it is nummerically more robust and has the same complexity. – Kaveh Vahedipour Apr 05 '19 at 09:40
`dsyevr` seems to work with larger matrices (n > 32766). I'm still waiting on the result ... – Wollmich Apr 05 '19 at 13:15

score 0 · Accepted Answer · answered Apr 18 '19 at 14:46

This seems to be a bug of LAPACKE_dsyevd. Using LAPACKE_dsyevr is working as well with larger matrices.

I've added the following lines to the MKL class:

    public static double[,] SymmetricEigRelativelyRobustRepresentations(double[,] a, out double[] w)
    {
        int n1 = a.GetLength(0);
        int n2 = a.GetLength(1);
        if (n1 != n2) throw new System.Exception("Matrix must be square");
        double[,] b = Copy(a);
        int matrix_layout = 101; // row-major arrays
        char jobz = 'V'; // eigenvalues and eigenvectors are computed
        char range = 'A'; // the routine computes all eigenvalues
        char uplo = 'U'; // a stores the upper triangular part of A 
        int n = n2;
        int lda = n;
        int vl = 0;
        int vu = 0;
        int il = 0;
        int iu = 0;
        double abstol = 0;
        int m = n;
        w = new double[n];
        double[,] z = new double[n, n];
        int ldz = n;
        int[] isuppz = new int[2 * n];
        int info = _mkl.LAPACKE_dsyevr(matrix_layout, jobz, range, uplo, n, b, lda, vl, vu, il, iu, abstol, ref m, w, z, ldz, isuppz);
        return z;
    }

and the following lines to the _mkl class:

    [DllImport(DLLName, CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)]
    internal static extern lapack_int LAPACKE_dsyevr(
        int matrix_layout, char jobz, char range, char uplo, lapack_int n, [In, Out] double[,] a, lapack_int lda,
        double vl, double vu, lapack_int il, lapack_int iu, double abstol, [In, Out] ref lapack_int m, [In, Out] double[] w,
        [In, Out] double[,] z, lapack_int ldz, [In, Out] lapack_int[] isuppz);

Intel MKL LAPACKE_dsyevd with n > 32766 --> Not enough memory to allocate work array in LAPACKE_dsyevd

2 Answers2