title: Why is matrix multiplication faster with numpy than with ctypes in Python?
body: <p>I was trying to figure out the fastest way to do matrix multiplication and tried 3 different ways:</p>
  
  <ul>
  <li>Pure python implementation: no surprises here.</li>
  <li>Numpy implementation using <code>numpy.dot(a, b)</code></li>
  <li>Interfacing with C using <code>ctypes</code> module in Python.</li>
  </ul>
  
  <p>This is the C code that is transformed into a shared library:</p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  
  void matmult(float* a, float* b, float* c, int n) {
      int i = 0;
      int j = 0;
      int k = 0;
  
      /*float* c = malloc(nay * sizeof(float));*/
  
      for (i = 0; i &lt; n; i++) {
          for (j = 0; j &lt; n; j++) {
              int sub = 0;
              for (k = 0; k &lt; n; k++) {
                  sub = sub + a[i * n + k] * b[k * n + j];
              }
              c[i * n + j] = sub;
          }
      }
      return ;
  }
  </code></pre>
  
  <p>And the Python code that calls it:</p>
  
  <pre><code>def C_mat_mult(a, b):
      libmatmult = ctypes.CDLL("./matmult.so")
  
      dima = len(a) * len(a)
      dimb = len(b) * len(b)
  
      array_a = ctypes.c_float * dima
      array_b = ctypes.c_float * dimb
      array_c = ctypes.c_float * dima
  
      suma = array_a()
      sumb = array_b()
      sumc = array_c()
  
      inda = 0
      for i in range(0, len(a)):
          for j in range(0, len(a[i])):
              suma[inda] = a[i][j]
              inda = inda + 1
          indb = 0
      for i in range(0, len(b)):
          for j in range(0, len(b[i])):
              sumb[indb] = b[i][j]
              indb = indb + 1
  
      libmatmult.matmult(ctypes.byref(suma), ctypes.byref(sumb), ctypes.byref(sumc), 2);
  
      res = numpy.zeros([len(a), len(a)])
      indc = 0
      for i in range(0, len(sumc)):
          res[indc][i % len(a)] = sumc[i]
          if i % len(a) == len(a) - 1:
              indc = indc + 1
  
      return res
  </code></pre>
  
  <p>I would have bet that the version using C would have been faster ... and I'd have lost ! Below is my benchmark which seems to show that I either did it incorrectly, or that <code>numpy</code> is stupidly fast:</p>
  
  <p><img src="https://i.stack.imgur.com/TLhDi.png" alt="benchmark"></p>
  
  <p>I'd like to understand why the <code>numpy</code> version is faster than the <code>ctypes</code> version, I'm not even talking about the pure Python implementation since it is kind of obvious.</p>