在尝试并行计算多个矩阵的特征值和特征向量时，我发现LAPACK的dsyevr函数似乎不是线程安全的。
- 有人知道这个问题吗？
- 我的代码有什么问题吗？（请参见下面的最小示例）
- 欢迎提供适用于密集矩阵的特征求解实现，该实现既不太慢，又绝对是线程安全的。
以下是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
The code was compiled with gcc 4.4.5 and linked against openblas (containing LAPACK) (libopenblas_sandybridge-r0.2.8.so). Another version of LAPACK was also tested. The program was tested by calling LAPACK directly from C (without LAPACKE), and the results were the same. Substituting
dsyevr with the dsyevd function (and adjusting arguments) had no effect. Here are the compilation instructions used: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 平台上完成的。该程序存在一致性问题。为了使内容更加通俗易懂，本文不进行过多解释，但会保留 HTML 标签。
icpc抱怨将void指针分配给另一种类型的指针，因此我不得不在solve函数中的malloc中添加(int*)等。 数值误差仅在并行部分中发生。 顺序部分对所有结果都恰好为零。 - Stefanicc（英特尔的C编译器）抱怨你不能在循环中定义变量（for (int i = 0...）。这就是我认为你在用C++的原因。 - Stefan