numpy, звонящий sse2 через ctypes

Короче говоря, я пытаюсь звонить в общую библиотеку из Python, более конкретно, от numpy. Общая библиотека реализована в C, использующем sse2 инструкции. Включая оптимизацию, т.е. создавая библиотеку с-O2 или –O1, я сталкиваюсь со странным segfaults при вызове в общую библиотеку через ctypes. Отключая оптимизацию (-O0), все удается как ожидалось, как имеет место при соединении библиотеки с c-программой непосредственно (оптимизированный или не). Присоединенный Вы находите отрезанный, который показывает очерченное поведение в моей системе. С включенной оптимизацией gdb сообщает о segfault в __ builtin_ia32_loadupd (__ P) в emmintrin.h:113. О значении __ P сообщают, как оптимизировано.

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

Флаги компилятора: gcc-o libtest.so - совместно использованный-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)

Вызов теста от C работает как ожидалось:

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;

}

Компиляция и вызов:
gcc-std=c99-otestc-msse2-L.-ltest call_from_c.c
экспортируйте $ LD_LIBRARY_PATH= {LD_LIBRARY_PATH}:.
./testc
... работы.

Моя система:

• 2.6.31-22-универсальный Ubuntu Linux i686
• Компилятор: gcc (4.4.1-4ubuntu9 Ubuntu)
• Python: Python 2.6.4 (r264:75706, 7 декабря 2009, 18:45:15) [GCC 4.4.1]
• Numpy: 1.4.0

Я взял условия (cf. код Python), что y выровненный, и выравнивание x не должно иметь значения (я думаю; явно выравнивание x не решает проблему хотя).

Обратите внимание также, что я использую _mm_loadu_pd вместо _mm_load_pd при загрузке b и f. Для версии C-only _mm_load_pd работы (как ожидалось). Однако при вызывании функции через ctypes, использующий _mm_load_pd всегда segfaults (независимый от оптимизации).

Я попробовал несколько дней для разбираний в этой проблеме без успеха..., и я нахожусь на грани, избивающей мой монитор до смерти. Любое входное приветствие. Daniel