Durante il tentativo di calcolare autovalori e autovettori di diversi matrici in parallelo, ho trovato che LAPACKs dsyevr funzione non sembra essere thread safe.Non dovrebbe LAPACKs funzione dsyevr (per autovalori e autovettori) essere thread-safe?
- È noto a tutti?
- c'è qualcosa di sbagliato con il mio codice? (Vedi esempio sotto il minimo)
- Eventuali suggerimenti di un'implementazione eigensolver per le matrici dense che non è troppo lento ed è sicuramente thread-safe è il benvenuto.
Ecco un esempio di codice minimo in C illustra il problema:
#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;
}
e questo è ciò che si ottiene (vedi differenza di blocco d'uscita scorso, che ti dice, che gli autovettori sono sbagliate, anche se autovalori sono ok):
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
Il codice è stato compilato con gcc 4.4.5 e collegati contro openblas (contenenti LAPACK) (libopenblas_sandybridge-r0.2.8.so), ma è stato anche testato con un'altra versione LAPACK. È stato anche testato il LAPACK direttamente da C (senza LAPACKE), stesso risultato. Sostituire dsyevr con la funzione dsyevd (e regolare gli argomenti) non ha avuto alcun effetto.
Infine, ecco le istruzioni di compilazione che ho usato:
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
Purtroppo Google non ha risposto alle mie domande, in modo che qualsiasi suggerimento è il benvenuto!
EDIT: Per assicurarsi che tutto è ok con le versioni BLAS e LAPACK ho preso il riferimento LAPACK (compresi BLAS e LAPACKE) da http://www.netlib.org/lapack/ (versione 3.4.2) Compilare il codice di esempio è stato un po ' difficile, ma ha finalmente lavorare con la compilazione separata e il collegamento:
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
la build del netLib LAPACK/BLAS e il programma di esempio è stato fatto su un Darwin 12.4.0 x86_64 e una piattaforma Linux 3.2.0-0.bpo.4-amd64 x86_64. È possibile osservare un comportamento errato coerente del programma.
La prima cosa da verificare è assicurarsi che si stiano veramente passando i set di dati disgiunti alle chiamate parallele. Secondo, LAPACK è implementato in FORTRAN, quindi dovresti in qualche modo ottenere quella lingua taggata nella domanda. – jxh
Buon consiglio.Ho modificato il codice di esempio per essere più leggibile e più sicuro con l'aritmetica del puntatore w.r.t. A tal fine, sono stati inclusi typedef per il vettore e la matrice. Ho anche controllato che i blocchi di memoria passati a LAPACK venissero disgiunti da alcune asserzioni aggiuntive (non incluse nel codice di esempio modificato per mantenerlo il più breve possibile). – dastrobu
http://www.netlib.org/lapack/lapack-3.3.0.html afferma che tutte le routine in LAPACK 3.3 sono thread-safe. –