複数の行列の固有値と固有ベクトルを並列に計算しようとしたときに、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 を直接呼び出すこともテストされ、同じ結果が得られました。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 は私の質問に答えてくれませんでした。ヒントがあれば大歓迎です!
編集: 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プラットフォーム上で行われました。プログラムの一貫した誤動作が観察されます。