Podczas próby obliczenia wartości własnych i wektorów własnych z kilku macierzy równolegle, stwierdziłem, że funkcja dsaevr LAPACK-ów nie wydaje się być bezpieczna dla wątków.Czy funkcja LAPACK dsyevr (dla wartości własnych i wektorów własnych) nie powinna być wątkowo bezpieczna?
- Czy to jest znane?
- Czy coś jest nie tak z moim kodem? (zobacz minimalny przykład poniżej)
- Wszelkie sugestie dotyczące implementacji eigensolver dla gęstych macierzy, które nie są zbyt wolne i są zdecydowanie bezpieczne dla wątków, są mile widziane.
Oto minimalne przykładowy kod w C, który demonstruje problem:
#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;
}
i to jest to, co masz (zobacz różnicę w ostatnim bloku wyjściowego, który mówi ci, że wektory własne są złe, jakkolwiek wartości własne są 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
kod został opracowany z gcc 4.4.5 i połączone na openblas (zawierających Lapack) (libopenblas_sandybridge-r0.2.8.so), lecz również testy z inną wersją Lapack. Wywołano także wywoływanie LAPACK bezpośrednio z C (bez LAPACKE), te same wyniki. Zastępowanie funkcji dsyevr przez funkcję dsyevd (i dostosowywanie argumentów) również nie przyniosło efektu.
Wreszcie, tutaj jest instrukcja kompilacji użyłem:
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
Niestety Google nie odpowiedział na moje pytania, więc każda wskazówka jest mile widziany!
EDIT: Aby upewnić się, że wszystko jest w porządku z wersjami Blas i LAPACK Wziąłem odniesienie LAPACK (łącznie BLAS i LAPACKE) z http://www.netlib.org/lapack/ (wersja 3.4.2) kompilacji kodu przykładem był trochę trudne, ale w końcu pracować z oddzielną kompilację i konsolidację:
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
kompilacji z netlib LAPACK/BLAS i przykładowego programu przeprowadzono na Darwin 12.4.0 x86_64 i platforma Linux 3.2.0-0.bpo.4-amd64 x86_64. Można zaobserwować konsekwentne złe zachowanie programu.
Pierwszą rzeczą, którą należy sprawdzić, jest upewnienie się, że rzeczywiście przekazywane są rozłączne zestawy danych do połączeń równoległych. Po drugie, LAPACK jest zaimplementowany w FORTRAN, więc powinieneś jakoś zaimportować ten język do pytania. – jxh
Dobra rada.Zmodyfikowałem przykładowy kod, aby był bardziej czytelny i bezpieczniejszy w arytmetyce wskaźnika w.r.t. W tym celu uwzględniono typedefs dla wektora i matrycy. Sprawdziłem również, że fragmenty pamięci przekazane do LAPACK-a są rozłączane przez dodatkowe asercje (nie zawarte w edytowanym kodzie przykładowym, aby zachować go jak najkrótszym). – dastrobu
http://www.netlib.org/lapack/lapack-3.3.0.html stwierdza, że wszystkie procedury w LAPACK 3.3 są bezpieczne dla wątków. –