title: Allocate intermediate multidimensional arrays in Cython without acquiring the GIL
body: <p>I'm trying to use Cython to parallelize an expensive operation which involves generating intermediate multidimensional arrays.</p>
  
  <p>The following very simplified code illustrates the sort of thing I'm trying to do:</p>
  
  <pre><code>import numpy as np
  cimport cython
  cimport numpy as np
  from cython.parallel cimport prange
  from libc.stdlib cimport malloc, free
  
  
  @cython.boundscheck(False)
  @cython.wraparound(False)
  def embarrasingly_parallel_example(char[:, :] A):
  
      cdef unsigned int m = A.shape[0]
      cdef unsigned int n = A.shape[1]
      cdef np.ndarray[np.float64_t, ndim = 2] out = np.empty((m, m), np.float64)
      cdef unsigned int ii, jj
      cdef double[:, :] tmp
  
      for ii in prange(m, nogil=True):
          for jj in range(m):
  
              # allocate a temporary array to hold the result of
              # expensive_function_1
              tmp_carray = &lt;double * &gt; malloc((n ** 2) * sizeof(double))
  
              # a 2D typed memoryview onto tmp_carray
              tmp = &lt;double[:n, :n] &gt; tmp_carray
  
              # shove the intermediate result in tmp
              expensive_function_1(A[ii, :], A[jj, :], tmp)
  
              # get the final (scalar) output for this ii, jj
              out[ii, jj] = expensive_function_2(tmp)
  
              # free the intermediate array
              free(tmp_carray)
  
      return out
  
  
  # some silly examples - the actual operation I'm performing is a lot more
  # involved
  # ------------------------------------------------------------------------
  @cython.boundscheck(False)
  @cython.wraparound(False)
  cdef void expensive_function_1(char[:] x, char[:] y, double[:, :] tmp):
  
      cdef unsigned int m = tmp.shape[0]
      cdef unsigned int n = x.shape[0]
      cdef unsigned int ii, jj
  
      for ii in range(m):
          for jj in range(m):
              tmp[ii, jj] = 0
              for kk in range(n):
                  tmp[ii, jj] += (x[kk] + y[kk]) * (ii - jj)
  
  
  @cython.boundscheck(False)
  @cython.wraparound(False)
  cdef double expensive_function_2(double[:, :] tmp):
  
      cdef unsigned int m = tmp.shape[0]
      cdef unsigned int ii, jj
      cdef double result = 0
  
      for ii in range(m):
          for jj in range(m):
              result += tmp[ii, jj]
  
      return result
  </code></pre>
  
  <p>There seems to be at least two reasons why this fails to compile:</p>
  
  <ol>
  <li><p>Based on the output of <code>cython -a</code>, the creation of the typed memory view here:</p>
  
  <pre><code>cdef double[:, :] tmp = &lt;double[:n, :n] &gt; tmp_carray
  </code></pre>
  
  <p>seems to involve Python API calls, and I therefore can't release the GIL to allow the outer loop to run in parallel.</p>
  
  <p>I was under the impression that typed memory views were not Python objects, and therefore a child process ought to be able to create them without first acquiring the GIL. Is this the case?</p></li>
  </ol>
  
  <p><s>2. Even if I replace <code>prange(m, nogil=True)</code> with a normal <code>range(m)</code>, Cython still doesn't seem to like the presence of a <code>cdef</code> within the inner loop:</s></p>
  
  <pre><code>    Error compiling Cython file:
      ------------------------------------------------------------
      ...
                  # allocate a temporary array to hold the result of
                  # expensive_function_1
                  tmp_carray = &lt;double*&gt; malloc((n ** 2) * sizeof(double))
  
                  # a 2D typed memoryview onto tmp_carray
                  cdef double[:, :] tmp = &lt;double[:n, :n]&gt; tmp_carray
                      ^
      ------------------------------------------------------------
  
      parallel_allocate.pyx:26:17: cdef statement not allowed here
  </code></pre>
  
  <h2>Update</h2>
  
  <p>It turns out that the second problem was easily solved by moving</p>
  
  <pre><code> cdef double[:, :] tmp
  </code></pre>
  
  <p>outside of the <code>for</code> loop, and just assigning</p>
  
  <pre><code> tmp = &lt;double[:n, :n] &gt; tmp_carray
  </code></pre>
  
  <p>within the loop. I still don't fully understand why this is necessary, though.</p>
  
  <p>Now if I try to use <code>prange</code> I hit the following compilation error:</p>
  
  <pre><code>Error compiling Cython file:
  ------------------------------------------------------------
  ...
              # allocate a temporary array to hold the result of
              # expensive_function_1
              tmp_carray = &lt;double*&gt; malloc((n ** 2) * sizeof(double))
  
              # a 2D typed memoryview onto tmp_carray
              tmp = &lt;double[:n, :n]&gt; tmp_carray
                 ^
  ------------------------------------------------------------
  
  parallel_allocate.pyx:28:16: Memoryview slices can only be shared in parallel sections
  </code></pre>