numpy appeler sse2 via ctypes

Question

En bref, j'essaie d'appeler une librairie partagée à partir de python, plus précisément de numpy. La bibliothèque partagée est implémentée en C en utilisant les instructions sse2. En activant l'optimisation, c'est-à-dire la construction de la bibliothèque avec -O2 ou -O1, je suis confronté à d'étranges segfaults lorsque j'appelle dans la bibliothèque partagée via des ctypes. Désactivation de l'optimisation (-O0), tout se passe comme prévu, comme c'est le cas lors de la liaison directe de la bibliothèque à un c-programme (optimisé ou non). Vous trouverez ci-joint un extrait qui présente le comportement délimité sur mon système. Avec l'optimisation activée, gdb signale un segfault dans __builtin_ia32_loadupd (__P) à emmintrin.h: 113. La valeur de __P est signalée comme optimisée.

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

drapeaux du compilateur: 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) 

appel test à partir de C fonctionne comme prévu:

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; 

} 

Compile et appelez:
gcc std = c99 -otestc -msse2 -L. -ltest call_from_c.c
export LD_LIBRARY_PATH = $ {LD_LIBRARY_PATH} :.
./testc
... fonctionne.

Mon système:

• Ubuntu Linux i686 2.6.31-22-générique
• compilateur: gcc (Ubuntu 4.4.1-4ubuntu9)
• Python: Python 2.6.4 (R264: 75706 7 Déc 2009 18:45:15) [GCC 4.4.1]
• Numpy: 1.4.0

j'ai prises provisions (cf. code python) qui y est aligné et l'alignement de x ne devrait pas mat r (je pense; aligner x explicitement ne résout pas le problème). Notez également que j'utilise _mm_loadu_pd au lieu de _mm_load_pd lors du chargement de b et f. Pour la version C uniquement, _mm_load_pd fonctionne (comme prévu). Toutefois, lors de l'appel de la fonction via les types de code à l'aide de _mm_load_pd toujours segfaults (indépendant de l'optimisation).

J'ai essayé plusieurs jours de résoudre ce problème sans succès ... et je suis sur le point de battre mon écran à mort. Toute entrée bienvenue. Daniel

Answer 1

Essayez de mettre à niveau vers Numpy 1.5.0b2. Il suffit d'exécuter ce qui suit (mais attention il peut briser d'autres choses (vous devrez recompiler tous pyrex):

sudo easy_install -U numpy 

J'avais des problèmes similaires avec ctypes quand je tentais d'utiliser H5PY (je devais recompiler le .deb pour obtenir avec la dernière version de numpy) et aussi il y avait des problèmes majeurs avec weave que la dernière mise à jour corrigée.

Answer 2

Essayez de construire votre extension en utilisant le système de construction numpy à rabais cflags potentiels/ldflags différences: http://projects.scipy.org/numpy/wiki/NumpySconsExtExamples

Answer 3

Je viens de me mordu par ce essayant d'appeler un certain ESS code de python, le problème semble être que GCC veut pour supposer que la pile est alignée sur des limites de 16 octets (le type natif le plus grand sur l'architecture, c'est-à-dire les types SSE), et calcule tous les décalages avec cette hypothèse. Lorsque cette hypothèse est fausse, les instructions SSE seront piégées.

La réponse semble être de compiler avec

gcc -mstackrealign

qui change la fonction prologues pour toujours aligner la pile à 16 octets.