【发布时间】：2013-08-13 18:17:31
【问题描述】：
在尝试并行计算 几个 矩阵的特征值和特征向量时，我发现 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
代码使用 gcc 4.4.5 编译并与 openblas（包含 LAPACK）（libopenblas_sandybridge-r0.2.8.so）链接，但也使用另一个 LAPACK 版本进行了测试。还测试了直接从 C（没有 LAPACKE）调用 LAPACK，结果相同。用dsyevd 函数替换dsyevr（并调整参数）也没有效果。
最后，这是我使用的编译指令：
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 平台上完成的。可以观察到程序的一贯不当行为。
【问题讨论】：
-
首先要检查的是确保您确实将不相交的数据集传递给并行调用。其次，LAPACK 是在 FORTRAN 中实现的，因此您应该以某种方式将该语言标记到问题中。
-
好建议。我修改了示例代码，使其更具可读性和更安全的 w.r.t 指针算法。为此，包括了向量和矩阵的类型定义。我还检查了传递给 LAPACK 的内存块是否被一些额外的断言分开（未包含在编辑的示例代码中以使其尽可能短）。
-
netlib.org/lapack/lapack-3.3.0.html 声明 LAPACK 3.3 中的所有例程都是线程安全的。
-
我用英特尔 MKL 和 icpc 尝试了您的代码，在第二次运行时，我只收到数字错误 (malloc 语句并将头文件更改为
mkl_lapacke.h... -
“必须添加类型...”我的意思是，
icpc抱怨将void指针分配给另一种类型的指针。所以我不得不添加例如(int*)到solve中的malloc。数值误差仅发生在平行部分。对于所有结果，顺序部分正好为零。
标签： c thread-safety fortran lapack eigenvector