numpy calling sse2 via ctypes

Question

In breve, sto provando a chiamare in una libreria condivisa da python, in particolare da numpy. La libreria condivisa è implementata in C usando le istruzioni sse2. Abilitando l'ottimizzazione, cioè costruendo la libreria con -O2 o -O1, sto affrontando strani segoults quando si chiama nella libreria condivisa tramite ctypes. Disattivando l'ottimizzazione (-O0), tutto funziona come previsto, come nel caso del collegamento diretto della libreria a un programma c (ottimizzato o meno). In allegato trovi un taglio che mostra il comportamento delineato sul mio sistema. Con l'ottimizzazione abilitata, gdb riporta un segfault in __builtin_ia32_loadupd (__P) su emmintrin.h: 113. Il valore di __P è segnalato come ottimizzato.

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; 
} 

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

test. 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) 

Chiamata di prova dalla C funziona come previsto:

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; 

} 

Compilare e chiamare:
gcc -std = c99 -otestc -msse2 -L. -ltest call_from_c.c
export LD_LIBRARY_PATH = $ {LD_LIBRARY_PATH} :.
./testc
... funziona.

Il mio sistema:

• Ubuntu Linux i686 2.6.31-22-generic
• compilatore gcc (Ubuntu 4.4.1-4ubuntu9)
• Python: Python 2.6.4 (r264: 75706 7 Dicembre 2009 18:45:15) [GCC 4.4.1]
• Numpy: disposizioni 1.4.0

ho preso (cfr codice pitone) che y è allineato e l'allineamento di x non dovrebbe opacizzare r (penso; allineare esplicitamente x non risolve il problema però).

Nota anche che io uso _mm_loadu_pd invece di _mm_load_pd quando carico b e f. Per la versione solo C _mm_load_pd funziona (come previsto). Tuttavia, quando si chiama la funzione tramite ctypes usando _mm_load_pd sempre segfaults (indipendentemente dall'ottimizzazione).

Ho provato diversi giorni per risolvere questo problema senza successo ... e sono sul punto di battere il mio monitor a morte. Qualsiasi input benvenuto. Daniel

Answer 1

Hai provato a eseguire l'aggiornamento a Numpy 1.5.0b2. Basta eseguire il seguente (ma attenzione che potrebbe rompere le altre cose (si dovrà ricompilare tutto pyrex):

sudo easy_install -U numpy 

ho avuto problemi simili con ctypes quando stavo cercando di usare H5PY (ho dovuto ricompilare il .deb per ottenere l'ultima versione di numpy) e inoltre ci sono stati grossi problemi con weave che l'aggiornamento più recente è stato corretto.

Answer 2

Provare a costruire la propria estensione utilizzando sistema di compilazione NumPy a scontare i potenziali CFLAGS/LDFLAGS differenze: http://projects.scipy.org/numpy/wiki/NumpySconsExtExamples

Answer 3

Ho appena ricevuto morso da questo cercando di chiamare un certo SSE-codice python, il problema sembra essere che GCC vuole per assumere che lo stack sia allineato su limiti di 16 byte (il più grande tipo nativo sull'architettura, cioè i tipi SSE), e calcola tutto l'offset con tale ipotesi. Quando questa ipotesi è falsa, le istruzioni SSE si intrappolano.

La risposta sembra essere quella di compilare con

gcc -mstackrealign

che modifica i prologhi delle funzioni per allineare sempre lo stack a 16 byte.