¿La función dsyev LAPACKs dsyevr(para valores propios y eigenvectores) no debe ser segura para subprocesos?

thread-safety fortran (4)

Al tratar de calcular valores propios y vectores propios de varias matrices en paralelo, encontré que la función dsyev LAPACKs no parece ser segura para subprocesos.

¿Esto es conocido por alguien?
¿Hay algún problema con mi código? (vea el ejemplo mínimo a continuación)
Cualquier sugerencia de una implementación de eigensolver para matrices densas que no sea demasiado lenta y definitivamente sea segura para subprocesos es bienvenida.

Aquí hay un ejemplo de código mínimo en C que demuestra el 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; }

y esto es lo que obtienes (mira la diferencia en el último bloque de salida, que te dice que los autovectores están equivocados, aunque los autovalores están bien):

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

El código se compiló con gcc 4.4.5 y se vinculó con openblas (que contiene LAPACK) (libopenblas_sandybridge-r0.2.8.so) pero también se probó con otra versión de LAPACK. Llamar a LAPACK directamente desde C (sin LAPACKE) también se probó, los mismos resultados. Sustituir dsyevr por la función dsyevd (y ajustar argumentos) tampoco tuvo ningún efecto.

Finalmente, aquí está la instrucción de compilación que utilicé:

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

Lamentablemente, Google no respondió mis preguntas, por lo que cualquier sugerencia es bienvenida.

EDITAR: para asegurarse de que todo está bien con las versiones BLAS y LAPACK tomé el LAPACK de referencia (incluyendo BLAS y LAPACKE) desde http://www.netlib.org/lapack/ (versión 3.4.2) Compilando el código de ejemplo fue un poco complicado, pero finalmente funcionó con la compilación y el enlace por separado:

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 compilación de netlib LAPACK / BLAS y el programa de ejemplo se realizaron en una plataforma Darwin 12.4.0 x86_64 y Linux 3.2.0-0.bpo.4-amd64 x86_64 . Se puede observar una mala conducta consistente del programa.

Finalmente recibí la explicación del equipo LAPACK, que me gustaría citar (con permiso):

Creo que el problema que está viendo puede deberse a la compilación de la versión de FORTRAN de la biblioteca LAPACK que está utilizando. Utilizando gfortran 4.8.0 (en Mac OS 10.8), podría reproducir el problema que vi si compilo LAPACK con la opción -O3 para gfortran. Si recompilo el LAPACK y hago referencia a la biblioteca BLAS con -fopenmp -O3, el problema desaparece. Hay una nota en el manual de gfortran que dice "-fopenmp implica -frecursive, es decir, todas las matrices locales se asignarán en la pila", por lo que puede haber matrices locales usadas en algunas rutinas auxiliares llamadas por dsyevr para las cuales la configuración predeterminada de la el compilador hace que se asignen de una manera que no sea segura para subprocesos. En cualquier caso, la asignación de estos en la pila, lo que parece hacer -fopenmp, abordará este problema.

Puedo confirmar que esto resuelve el problema para netlib-BLAS / LAPACK. Se debe tener en cuenta que el tamaño de la pila es limitado y posiblemente se debe ajustar si las matrices se agrandan y / o muchas.

OpenBLAS debe compilarse con USE_OPENMP=1 y USE_THREAD=1 para obtener una biblioteca con un único subproceso y una seguridad de subprocesos.

Con estos indicadores de compilación y marca, el programa de ejemplo se ejecuta correctamente con todas las bibliotecas. Sigue siendo una pregunta abierta, cómo uno se asegura de que el usuario a quien uno le entrega el código al final esté vinculándose a una biblioteca BLAS / LAPACK correctamente compilada. Si el usuario acaba de obtener una falla de segmentación, se podría agregar una nota en un archivo README, pero como el error es más sutil, ni siquiera se garantiza que el usuario reconozca el error (los usuarios no leen el archivo README por defecto ;-) ). No es una buena idea enviar un BLAS / LAPACK con el código de uno, ya que es la idea básica de BLAS / LAPACK que todos tengan una versión optimizada específicamente para su computadora. Las ideas son bienvenidas ...

Parece que les pides a los desarrolladores de LAPACK que introduzcan una "solución". De hecho, agregaron -frecursive a los indicadores del compilador en make.inc.example. Al probar su código de ejemplo, la corrección parece irrelevante (los errores numéricos no desaparecen) y no es deseable (un posible golpe de rendimiento).

Incluso si la corrección era correcta, -frecursive está implícito en -fopenmp, por lo que las personas que usan indicadores consistentes están en el lado seguro (los que usan software preempaquetado no lo son).

Para concluir, corrija su código en lugar de confundir a las personas.

Re otra biblioteca: GSL. Es threadsafe Pero eso significa que debe crear espacios de trabajo para cada subproceso y asegurarse de que cada subproceso lo utilice en el espacio de trabajo, por ejemplo, punteros de índice por número de subproceso.

[la siguiente respuesta se agregó antes de que se conociera la solución correcta]

Descargo de responsabilidad: La siguiente es una solución, ni resuelve el problema original, ni explica qué pasa con LAPACK. Sin embargo, puede ayudar a las personas que enfrentan el mismo problema.

La vieja versión f2c de LAPACK, llamada CLAPACK, no parece tener el problema de seguridad de hilos. Tenga en cuenta que esta no es una interfaz C para la biblioteca Fortran, sino una versión de LAPACK que se ha traducido automáticamente a C.

Compilando y vinculándolo con la última versión de CLAPACK (3.2.1) funcionó. Entonces CLAPACK parece ser seguro para hilos a este respecto. Por supuesto, el rendimiento de CLAPACK no es el de netlib-BLAS / LAPACK o incluso el de OpenBLAS / LAPACK, pero al menos no es tan malo como el de GSL.

Aquí hay algunos tiempos para todas las variantes LAPACK probadas (y GSL) para la diagonalización de una matriz de 1000 x 1000 (en una secuencia de curso) inicializada con la función init (vea la pregunta para la definición).

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

Esto indica que GSL definitivamente no es una buena solución para matrices grandes con una dimensión de unos pocos miles, especialmente si tiene muchas de ellas.