Должна ли функция LAPACK dsyevr (для собственных значений и собственных векторов) быть потокобезопасной?

При попытке вычислить собственные значения и собственные векторы из нескольких матриц параллельно, я обнаружил, что функция dsevr LAPACKs не кажется потокобезопасной.

• Известно ли это кому-либо?
• Что-то не так с моим кодом? (см. минимальный пример ниже)
• Любые предложения реализации eigensolver для плотных матриц, которые не слишком медленны и, безусловно, являются потокобезопасными, приветствуются.

Вот пример минимального кода в C, который демонстрирует проблему:

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <omp.h>
#include "lapacke.h"

#define M 8 /* number of matrices to be diagonalized */
#define N 1000 /* size of each matrix (real, symmetric) */

typedef double vec_t[N]; /* type for length N vector */
typedef double mtx_t[N][N]; /* type for N x N matrices */

void 
init(int m, int n, mtx_t *A){
    /* init m symmetric n x x matrices */
    srand(0);
    for (int i = 0; i < m; ++i){
        for (int j = 0; j < n; ++j){
            for (int k = 0; k <= j; ++k){
                A[i][j][k] = A[i][k][j] = (rand()%100-50) / (double)100.;
            }
        }
    }
}

void 
solve(int n, double *A, double *E, double *Q){
    /* diagonalize one matrix */
    double tol = 0.;
    int *isuppz = malloc(2*n*sizeof(int)); assert(isuppz);
    int k;
    int info = LAPACKE_dsyevr(LAPACK_COL_MAJOR, 'V', 'A', 'L', 
                              n, A, n, 0., 0., 0, 0, tol, &k, E, Q, n, isuppz);
    assert(!info);
    free(isuppz);
}

void 
s_solve(int m, int n, mtx_t *A, vec_t *E, mtx_t *Q){
    /* serial solve */
    for (int i = 0; i < m; ++i){
        solve(n, (double *)A[i], (double *)E[i], (double *)Q[i]);
    }
}

void 
p_solve(int m, int n, mtx_t *A, vec_t *E, mtx_t *Q, int nt){
    /* parallel solve */
    int i;
    #pragma omp parallel for schedule(static) num_threads(nt) \
        private(i) \
        shared(m, n, A, E, Q)
    for (i = 0; i < m; ++i){
        solve(n, (double *)A[i], (double *)E[i], (double *)Q[i]);
    }
}

void 
analyze_results(int m, int n, vec_t *E0, vec_t *E1, mtx_t *Q0, mtx_t *Q1){
    /* compare eigenvalues */
    printf("\nmax. abs. diff. of eigenvalues:\n");
    for (int i = 0; i < m; ++i){
        double t, dE = 0.;
        for (int j = 0; j < n; ++j){
            t = fabs(E0[i][j] - E1[i][j]);
            if (t > dE) dE = t;
        }
        printf("%i: %5.1e\n", i, dE);
    }

    /* compare eigenvectors (ignoring sign) */
    printf("\nmax. abs. diff. of eigenvectors (ignoring sign):\n");
    for (int i = 0; i < m; ++i){
        double t, dQ = 0.;
        for (int j = 0; j < n; ++j){
            for (int k = 0; k < n; ++k){
                t = fabs(fabs(Q0[i][j][k]) - fabs(Q1[i][j][k]));
                if (t > dQ) dQ = t;
            }
        }
        printf("%i: %5.1e\n", i, dQ);
    }
}


int main(void){
    mtx_t *A = malloc(M*N*N*sizeof(double)); assert(A);
    init(M, N, A);

    /* allocate space for matrices, eigenvalues and eigenvectors */
    mtx_t *s_A = malloc(M*N*N*sizeof(double)); assert(s_A);
    vec_t *s_E = malloc(M*N*sizeof(double));   assert(s_E);
    mtx_t *s_Q = malloc(M*N*N*sizeof(double)); assert(s_Q);

    /* copy initial matrix */
    memcpy(s_A, A, M*N*N*sizeof(double));

    /* solve serial */
    s_solve(M, N, s_A, s_E, s_Q);

    /* allocate space for matrices, eigenvalues and eigenvectors */
    mtx_t *p_A = malloc(M*N*N*sizeof(double)); assert(p_A);
    vec_t *p_E = malloc(M*N*sizeof(double));   assert(p_E);
    mtx_t *p_Q = malloc(M*N*N*sizeof(double)); assert(p_Q);

    /* copy initial matrix */
    memcpy(p_A, A, M*N*N*sizeof(double));

    /* use one thread, to check that the algorithm is deterministic */
    p_solve(M, N, p_A, p_E, p_Q, 1); 

    analyze_results(M, N, s_E, p_E, s_Q, p_Q);

    /* copy initial matrix */
    memcpy(p_A, A, M*N*N*sizeof(double));

    /* use several threads, and see what happens */
    p_solve(M, N, p_A, p_E, p_Q, 4); 

    analyze_results(M, N, s_E, p_E, s_Q, p_Q);

    free(A);
    free(s_A);
    free(s_E);
    free(s_Q);
    free(p_A);
    free(p_E);
    free(p_Q);
    return 0;
}

и это то, что вы получаете (см. разницу в последнем блоке вывода, который говорит вам, что собственные векторы ошибочны, хотя собственные значения в порядке):

max. abs. diff. of eigenvalues:
0: 0.0e+00
1: 0.0e+00
2: 0.0e+00
3: 0.0e+00
4: 0.0e+00
5: 0.0e+00
6: 0.0e+00
7: 0.0e+00

max. abs. diff. of eigenvectors (ignoring sign):
0: 0.0e+00
1: 0.0e+00
2: 0.0e+00
3: 0.0e+00
4: 0.0e+00
5: 0.0e+00
6: 0.0e+00
7: 0.0e+00

max. abs. diff. of eigenvalues:
0: 0.0e+00
1: 0.0e+00
2: 0.0e+00
3: 0.0e+00
4: 0.0e+00
5: 0.0e+00
6: 0.0e+00
7: 0.0e+00

max. abs. diff. of eigenvectors (ignoring sign):
0: 0.0e+00
1: 1.2e-01
2: 1.6e-01
3: 1.4e-01
4: 2.3e-01
5: 1.8e-01
6: 2.6e-01
7: 2.6e-01

Код был скомпилирован с gcc 4.4.5 и связан с openblas (содержащий LAPACK) (libopenblas_sandybridge-r0.2.8.so), но также был протестирован с другой версией LAPACK. Вызов LAPACK непосредственно из C (без LAPACKE) также был протестирован с одинаковыми результатами. Подстановка dsyevr функцией dsyevd (и настройкой аргументов) также не имела эффекта.

Наконец, вот инструкция по компиляции, которую я использовал:

gcc -std=c99 -fopenmp -L/path/to/openblas/lib -Wl,-R/path/to/openblas/lib/ \
-lopenblas -lgomp -I/path/to/openblas/include main.c -o main

К сожалению, Google не ответил на мои вопросы, поэтому приветствуем любой привет!

EDIT: Чтобы убедиться, что все в порядке с версиями BLAS и LAPACK, я взял ссылку LAPACK (включая BLAS и LAPACKE) из http://www.netlib.org/lapack/ (версия 3.4.2) Компиляция кода примера была немного сложной, но, наконец, работала с отдельной компиляцией и компоновкой:

gcc -c -std=c99 -fopenmp -I../lapack-3.4.2/lapacke/include \
    netlib_dsyevr.c -o netlib_main.o
gfortran netlib_main.o ../lapack-3.4.2/liblapacke.a \
    ../lapack-3.4.2/liblapack.a ../lapack-3.4.2/librefblas.a \
    -lgomp -o netlib_main

Строка netlib LAPACK/BLAS и примерная программа были выполнены на платформе Darwin 12.4.0 x86_64 и Linux 3.2.0-0.bpo.4-amd64 x86_64. Можно наблюдать последовательное неправильное поведение программы.

Ответы

time -p ./gsl
real 17.45
user 17.42
sys 0.01

time -p ./netlib_dsyevr
real 0.67
user 0.84
sys 0.02

time -p ./openblas_dsyevr
real 0.66
user 0.46
sys 0.01

time -p ./clapack_dsyevr
real 1.47
user 1.46
sys 0.00