title: numpy calling sse2 via ctypes
body: <p>In brief, I am trying to call into a shared library from python, more specifically, from numpy. The shared library is implemented in C using sse2 instructions. Enabling optimisation, i.e. building the library with -O2 or –O1, I am facing strange segfaults when calling into the shared library via ctypes. Disabling optimisation (-O0), everything works out as expected, as is the case when linking the library to a c-program directly (optimised or not). Attached you find a snipped which exhibits the delineated behaviour on my system. With optimisation enabled, gdb reports a segfault in __builtin_ia32_loadupd (__P) at emmintrin.h:113. The value of __P is reported as optimised out.</p>
  
  <p>test.c:</p>
  
  <pre><code>#include &lt;emmintrin.h&gt;
  #include &lt;complex.h&gt;
  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 *) &amp;b );
  
      _p = (__m128d*) y;
  
      for(i=0; i&lt;m; ++i) {
          f = cexp(-I*x[i]);
          _f = _mm_loadu_pd( (double *) &amp;f );
          _x = _mm_loadu_pd( (double *) &amp;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;
  }
  </code></pre>
  
  <p>Compiler flags:
  gcc -o libtest.so -shared -std=c99 -msse2 -fPIC -O2 -g -lm test.c </p>
  
  <p>test.py:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>Calling test from C works as expected:</p>
  
  <p>call_from_c.c:</p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;complex.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;emmintrin.h&gt;
  
  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&lt;n; ++i) {
          x[i] = 1;
          y[i] = 0;
      }
  
      test(n, x, y);
      for(i=0; i&lt;n; ++i)
              printf("[%f %f]\n", creal(y[i]), cimag(y[i]));
  
      return 1;
  
  }
  </code></pre>
  
  <p>Compile and call:<br>
  gcc -std=c99 -otestc -msse2 -L. -ltest call_from_c.c<br>
  export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:.<br>
  ./testc<br>
  ... works.</p>
  
  <p>My system:</p>
  
  <ul>
  <li>Ubuntu Linux i686  2.6.31-22-generic </li>
  <li>Compiler: gcc (Ubuntu 4.4.1-4ubuntu9)</li>
  <li>Python: Python 2.6.4 (r264:75706, Dec 7 2009, 18:45:15) [GCC 4.4.1] </li>
  <li>Numpy: 1.4.0</li>
  </ul>
  
  <p>I have taken provisions (cf. python code) that y is aligned and the alignment of x should not matter (I think; explicitly aligning x does not solve the problem though).</p>
  
  <p>Note also that i use _mm_loadu_pd instead of _mm_load_pd when loading b and f. For the C-only version _mm_load_pd works (as expected). However, when calling the function via ctypes using
  _mm_load_pd always segfaults (independent of optimisation).</p>
  
  <p>I have tried several days to sort out this issue without success ... and I am on the verge beating my monitor to death. Any input welcome.
  Daniel</p>