llamadas numpy sse2 via ctypes

Question

En resumen, estoy tratando de llamar a una biblioteca compartida de python, más específicamente, de numpy. La biblioteca compartida se implementa en C usando las instrucciones sse2. Al habilitar la optimización, es decir, construir la biblioteca con -O2 o -O1, estoy enfrentando extraños errores al llamar a la biblioteca compartida a través de ctypes. Desactivando la optimización (-O0), todo funciona como se espera, como es el caso cuando se vincula la biblioteca a un programa c directamente (optimizado o no). Adjunto encuentra un recorte que exhibe el comportamiento delineado en mi sistema. Con la optimización habilitada, gdb informa una segfault en __builtin_ia32_loadupd (__P) en emmintrin.h: 113. El valor de __P se informa como optimizado.

test.c:

#include <emmintrin.h> 
#include <complex.h> 
void test(const int m, const double* x, double complex* y) { 

    int i; 
    __m128d _f, _x, _b; 
    double complex f __attribute__((aligned(16))); 
    double complex b __attribute__((aligned(16))); 
    __m128d* _p; 

    b = 1; 
    _b = _mm_loadu_pd((double *) &b); 

    _p = (__m128d*) y; 

    for(i=0; i<m; ++i) { 
     f = cexp(-I*x[i]); 
     _f = _mm_loadu_pd((double *) &f); 
     _x = _mm_loadu_pd((double *) &x[i]);  
     _f = _mm_shuffle_pd(_f, _f, 1); 
     *_p = _mm_add_pd(*_p, _f); 
     *_p = _mm_add_pd(*_p, _x); 
     *_p = _mm_mul_pd(*_p,_b); 
     _p++; 
    } 
    return; 
} 

banderas del compilador: gcc -o libtest.so -shared -std = c99 -msse2 -fPIC -O2 -g -lm test.c

prueba. PY:

import numpy as np 
import os 

def zerovec_aligned(nr, dtype=np.float64, boundary=16): 
    '''Create an aligned array of zeros. 
    ''' 
    size = nr * np.dtype(dtype).itemsize 
    tmp = np.zeros(size + boundary, dtype=np.uint8) 
    address = tmp.__array_interface__['data'][0] 
    offset = boundary - address % boundary 
    return tmp[offset:offset + size].view(dtype=dtype) 


lib = np.ctypeslib.load_library('libtest', '.') 
lib.test.restype = None 
lib.test.argtypes = [np.ctypeslib.ctypes.c_int, 
        np.ctypeslib.ndpointer(np.float64, flags=('C', 'A')), 
        np.ctypeslib.ndpointer(np.complex128, flags=('C', 'A', 'W'))] 


n = 13 
y = zerovec_aligned(n, dtype=np.complex128) 
x = np.ones(n, dtype=np.float64) 
# x = zerovec_aligned(n, dtype=np.float64) 
# x[:] = 1. 

lib.test(n,x,y) 

Llamada de prueba de C funciona como se esperaba:

call_from_c.c:

#include <stdio.h> 
#include <complex.h> 
#include <stdlib.h> 
#include <emmintrin.h> 

void test(const int m, const double* x, double complex* y); 

int main() { 

    int i; 
    const int n = 11; 
    double complex *y = (double complex*) _mm_malloc(n*sizeof(double complex), 16); 
    double *x = (double *) malloc(n*sizeof(double)); 
    for(i=0; i<n; ++i) { 
     x[i] = 1; 
     y[i] = 0; 
    } 

    test(n, x, y); 
    for(i=0; i<n; ++i) 
      printf("[%f %f]\n", creal(y[i]), cimag(y[i])); 

    return 1; 

} 

Compilar y llamar:
gcc -std = c99 -otestc -msse2 -L. -ltest call_from_c.c
export LD_LIBRARY_PATH = $ {LD_LIBRARY_PATH} :.
./testc
... funciona.

Mi sistema: i686

• Ubuntu Linux 2.6.31-22-generic
• Compilador: gcc (Ubuntu 4.4.1-4ubuntu9)
• Python: Python 2.6.4 (R264: 75706 7 Dic de 2009, 18:45:15) [GCC 4.4.1]
• Numpy: disposiciones 1.4.0

que he tomado (código Python cf.) que y es alineado y la alineación de x no debería ser mate r (creo; alinear explícitamente x no resuelve el problema, sin embargo).

Tenga en cuenta también que uso _mm_loadu_pd en lugar de _mm_load_pd al cargar b y f. Para la versión solo en C, _mm_load_pd funciona (como se esperaba). Sin embargo, al llamar a la función a través de ctypes utilizando _mm_load_pd siempre segfaults (independientemente de la optimización).

He intentado varios días resolver este problema sin éxito ... y estoy a punto de apalear mi monitor hasta la muerte. Cualquier entrada de bienvenida. Daniel

Answer 1

¿Has intentado actualizar a Numpy 1.5.0b2. Simplemente ejecute el siguiente (pero tenga cuidado de que podría romper otras cosas (que tendrá que volver a compilar todo pyrex):

sudo easy_install -U numpy 

que estaba teniendo problemas similares con ctypes cuando yo estaba tratando de usar H5PY (tuve que volver a compilar el .deb para obtener la última versión de numpy) y también hubo problemas importantes con el tejido que se corrigió en la última actualización.

Answer 2

Intenta construir su extensión usando el sistema de acumulación numpy descontar cflags potenciales/ldflags diferencias: http://projects.scipy.org/numpy/wiki/NumpySconsExtExamples

Answer 3

yo sólo le picaron por esta tratando de llamar a algunos SSE-código del pitón, el problema parece ser que GCC quiere suponer que la pila está alineada en los límites de 16 bytes (el tipo nativo más grande en la arquitectura, es decir, los tipos SSE), y calcula todo el desplazamiento con esa suposición. Cuando esa suposición es falsa, las instrucciones SSE se atraparán.

La respuesta parece ser la de compilar con

gcc -mstackrealign

que cambia los prólogos de función para alinear siempre la pila a 16 bytes.