title: Using memset instead of np.zeros in a Cython script for speed gains
body: <p>I started working on a SciPy interface to Fortran libraries (BLAS/LAPACK) as one can see here: <a href="https://stackoverflow.com/questions/44710838/calling-blas-lapack-directly-using-the-scipy-interface-and-cython">Calling BLAS / LAPACK directly using the SciPy interface and Cython</a> and came up with a solution but had to resort to using <code>numpy.zeros</code> which in effect killed any speed gains from calling the Fortran code directly.  The issue was the Fortran code needs an 0 valued output matrix (it operates on the matrix in-place in memory) to match the Numpy version (<code>np.outer</code>).  </p>
  
  <p>So I was a bit perplexed since a 1000x1000 zeros matrix in Python only take 8us (using %timeit, or 0.008ms) so why did adding Cython code kill the runtime, noting I'm also creating it on a memoryview? (basically it goes from 3ms to 8ms or so on a 1000 x 1000 matrix multiplication).  Then after searching around on SO I found a solution elsewhere using <code>memset</code> as the fastest array changing mechanism.  So I rewrote the whole code from the referenced post to the newer <code>memoryview</code> format and got to something like this (Cython <code>cyblas.pyx</code> file):</p>
  
  <pre><code>import cython
  import numpy as np
  cimport numpy as np
  from libc.string cimport memset #faster than np.zeros
  
  REAL = np.float64
  ctypedef np.float64_t REAL_t
  ctypedef np.uint64_t  INT_t
  
  cdef int ONE = 1
  cdef REAL_t ONEF = &lt;REAL_t&gt;1.0
  
  @cython.boundscheck(False)
  @cython.wraparound(False)
  @cython.cdivision(True)
  cpdef outer_prod(double[::1] _x, double[::1] _y, double[:, ::1] _output):
  
      cdef int M = _y.shape[0]
      cdef int N = _x.shape[0]
      memset(&amp;_output[0,0], 0, M*N)
  
      with nogil:
          dger(&amp;M, &amp;N, &amp;ONEF, &amp;_y[0], &amp;ONE, &amp;_x[0], &amp;ONE, &amp;_output[0,0], &amp;M)
  </code></pre>
  
  <h1>TEST SCRIPT</h1>
  
  <pre><code>import numpy as np;
  from cyblas import outer_prod;
  a=np.random.randint(0,100, 1000);
  b=np.random.randint(0,100, 1000);
  a=a.astype(np.float64)
  b=b.astype(np.float64)
  cy_outer=np.zeros((a.shape[0],b.shape[0]));
  np_outer=np.zeros((a.shape[0],b.shape[0]));
  
  %timeit outer_prod(a,b, cy_outer)
  %timeit np.outer(a,b, np_outer)
  </code></pre>
  
  <p>So this fixed my issue of resetting the output matrix values BUT ONLY TO ROW 125, BEYOND THAT THIS ISSUE STILL EXISTS (resolved see below).  I thought setting the <code>memset</code> length parameter to M*N would clear 1000*1000 in memory but apparently not.</p>
  
  <p>Does anyone have an idea on how to reset the whole output matrix to 0 using <code>memset</code>?  Much appreciated.</p>