不应该LAPACK dsyevr函数(对于特征值和特征向量)是线程安全的吗？

Question

在尝试并行计算几个矩阵的特征值和特征向量时,我发现LAPACK dsyevr函数似乎不是线程安全的.

• 这对任何人都知道吗？
• 我的代码有问题吗？(见下面的最小例子)
• 任何关于密集矩阵的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版本进行了测试.还测试了直接从C(没有LAPACKE)调用LAPACK,结果相同.代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

不幸的是谷歌没有回答我的问题,所以任何提示都是受欢迎的!

编辑: 为了确保BLAS和LAPACK版本的一切正常我从http://www.netlib.org/lapack/(版本3.4.2)参考LAPACK(包括BLAS和LAPACKE )编译示例代码有点棘手,但最终使用单独的编译和链接:

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平台上完成的.可以观察到程序的一致性错误行为.

Answer 1

我终于收到了LAPACK团队的解释,我想引用(经许可):

我认为您遇到的问题可能是由于您正在使用的LAPACK库的FORTRAN版本是如何编译的.使用gfortran 4.8.0(在Mac OS 10.8上),如果我用gfortran的-O3选项编译LAPACK,我可以重现你看到的问题.如果我重新编译LAPACK并使用-fopenmp -O3引用BLAS库,问题就会消失.gfortran手册中有一条说明"-fopenmp暗示-frecursive,即所有本地数组都将在堆栈上分配",因此在dsyevr调用的某些辅助例程中可能会使用本地数组,其默认设置为编译器使它们以非线程安全的方式分配.在任何情况下,在堆栈上分配这些 - -fopenmp似乎都会解决这个问题.

我可以确认这解决了netlib-BLAS/LAPACK的问题.应该记住,堆栈大小是有限的,并且如果矩阵变大和/或很多,可能需要调整.

必须使用USE_OPENMP=1和编译OpenBLAS USE_THREAD=1以获得单线程和线程安全的库.

使用这些编译器和make标志,示例程序可以正确运行所有库.这仍然是一个悬而未决的问题,如何确保最终一手牌代码的用户链接到正确编译的BLAS/LAPACK库？如果用户只是得到分段错误,可以在README文件中添加注释,但由于错误更加微妙,甚至不能保证用户可以识别错误(用户不通过以下方式读取自述文件)默认;-)).使用一些代码发送BLAS/LAPACK并不是一个好主意,因为BLAS/LAPACK的基本思想是每个人都有一个专门为他的计算机优化的版本.欢迎提出意见......