Intel MKL LAPACKE_dsyevd с n> 32766 -> Недостаточно памяти для выделения рабочего массива в LAPACKE_dsyevd

Question

Я хочу вычислить все собственные значения и все собственные векторы реальной симметричной матрицы, используя LAPACKE_dsyevd из Intel MKL (обновление 2019 года).

У меня следующий метод в 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;
    }
}

с

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);
}

и следующим тестовым кодом:

    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);

с

    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;
    }

Если n большечем 32766, происходит сбой со следующим сообщением об ошибке:

Недостаточно памяти для выделения рабочего массива в LAPACKE_dsyevd

Мой компьютер имеет 128 ГБ ОЗУ, чего должно быть достаточно.Мне известно о <gcAllowVeryLargeObjects enabled="true" />, и я установил его на true.Мне также известно о пределе 65535 ^ 2 для многомерного массива в C #, см. 2d-массив с более чем 65535 ^ 2 элементами -> Размеры массива превысили поддерживаемый диапазон .

Кстати, я могу вычислить разложение по собственным значениям, используя MATLAB для матриц с n = 40000 или больше.И я думаю, что MATLAB также использует Intel MKL в фоновом режиме для его вычисления.

Итак, как я могу вычислить разложение по собственным значениям очень больших матриц (n> 40000) в C # с использованием Intel MKL?

Answer 1

Кажется, это ошибка LAPACKE_dsyevd . Использование LAPACKE_dsyevr также работает с матрицами большего размера.

Я добавил следующие строки в класс MKL:

    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;
    }

и следующие строки для _mkl класса:

    [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);

Answer 2

Я не думаю, что это ваша проблема. Приведенное ниже определение предполагает, что w = new double[n]; в вашем случае слишком мало.

*  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.

Вы действительно должны всегда делать запрос рабочей области. Я знаю, что документация оставляет это открытым для пользователя, но это так удобно и помогает избежать вашей ситуации. Итак, вы знаете, как сделать запрос рабочей области? Если не ударить меня очень быстро.