Non dovrebbe LAPACKs funzione dsyevr (per autovalori e autovettori) essere thread-safe?

Question

Durante il tentativo di calcolare autovalori e autovettori di diversi matrici in parallelo, ho trovato che LAPACKs dsyevr funzione non sembra 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.

Answer 1

ho finalmente ricevuto la spiegazione da parte del team LAPACK, che vorrei citare (con autorizzazione):

Credo che il problema che state vedendo può essere causato da come la versione FORTRAN della libreria LAPACK stai usando è stato compilato. Usando gfortran 4.8.0 (su Mac OS 10.8), ho potuto riprodurre il problema che hai visto se compilavo LAPACK con l'opzione -O3 per gfortran. Se ricompongo il LAPACK e faccio riferimento alla libreria BLAS con -fopenmp -O3, il problema scompare. C'è una nota nel manuale di gfortran che dice "-fopenmp implica -frecursivo, cioè, tutti gli array locali saranno allocati nello stack", quindi potrebbero esserci matrici locali usate in alcune routine ausiliarie chiamate da dsyevr per le quali l'impostazione predefinita del il compilatore sta causando che vengano allocati in modo non thread-safe. In ogni caso, l'assegnazione di questi in pila, che sembra fare -fopenmp, risolverà questo problema.

Posso confermare che questo risolve il problema per netlib-BLAS/LAPACK. Si dovrebbe tenere presente che la dimensione dello stack è limitata e può essere corretta se le matrici diventano grandi e/o molte.

OpenBLAS deve essere compilato con USE_OPENMP=1 e USE_THREAD=1 per ottenere una libreria con thread singolo e thread-safe.

Con questi compilatori e contrassegni, il programma di esempio viene eseguito correttamente con tutte le librerie. Rimane una domanda aperta, in che modo ci si assicura che l'utente a cui appartiene il proprio codice utente alla fine si colleghi a una libreria BLAS/LAPACK compilata correttamente? Se l'utente ottiene un errore di segmentazione, è possibile aggiungere una nota in un file README, ma poiché l'errore è più sottile, non è nemmeno garantito che l'errore venga riconosciuto dall'utente (gli utenti non leggono il file README di default ;-)). Spedire un BLAS/LAPACK con un codice non è una buona idea, poiché è l'idea di base di BLAS/LAPACK che tutti abbiano una versione ottimizzata per il suo computer. Le idee sono benvenute ...

Answer 2

Re un'altra libreria: GSL. È infallibile. Ciò significa che devi creare spazi di lavoro per ciascun thread e assicurarti che ogni thread utilizzi lo spazio di lavoro, ad es. I puntatori dell'indice per numero di thread.

Answer 3

[la risposta che segue è stato aggiunto prima che la soluzione corretta è stata conosciuta]

responsabilità: La seguente è una soluzione, né si risolve il problema originale, né spiega che cosa va male con LAPACK. Può, tuttavia, aiutare le persone che affrontano lo stesso problema.

---

La vecchia versione f2c'ed di LAPACK, chiamato CLAPACK, non sembrano avere il problema di thread-sicurezza. Si noti che questa non è un'interfaccia C per la libreria fortran ma una versione di LAPACK che è stata automaticamente tradotta in C.

Compilare e collegarlo con l'ultima versione di CLAPACK (3.2.1) ha funzionato. Quindi CLAPACK sembra essere thread-safe in questo senso. Ovviamente, le prestazioni di CLAPACK non sono quelle di netlib-BLAS/LAPACK o anche di quelle di OpenBLAS/LAPACK, ma almeno non è così male come quella di GSL.

Ecco alcuni tempi per tutte le varianti LAPACK testate (e GSL) per la diagonalizzazione di una matrice 1000 x 1000 (su un thread ovviamente) inizializzata con la funzione init (vedere la domanda per definizione).

time -p ./gsl 
real 17.45 
user 17.42 
sys 0.01 

time -p ./netlib_dsyevr 
real 0.67 
user 0.84 
sys 0.02 

time -p ./openblas_dsyevr 
real 0.66 
user 0.46 
sys 0.01 

time -p ./clapack_dsyevr 
real 1.47 
user 1.46 
sys 0.00 

Questo indica che GSL è sicuramente non va bene soluzione per le grandi matrici di dimensione di alcune migliaia di persone, soprattutto se si dispone di molti di loro.

Answer 4

Sembra che tu abbia richiesto agli sviluppatori di LAPACK di introdurre una "correzione". Infatti, hanno aggiunto -comprensivo ai flag del compilatore in make.inc.example. Dalla verifica del codice di esempio, la correzione sembra irrilevante (gli errori numerici non scompaiono) e non è auspicabile (un possibile colpo alla performance).

Anche se la correzione era corretta, -frecursiva è implicita da -fopenmp, quindi le persone che usano flag coerenti sono sul lato sicuro (quelli che usano software preconfezionato non lo sono).

Per concludere, correggere il codice anziché confondere le persone.