Asigne matrices multidimensionales intermedias en Cython sin adquirir el GIL

Estoy tratando de usar Cython para paralelizar una operación costosa que implica generar matrices multidimensionales intermedias.

El siguiente código muy simplificado ilustra el tipo de cosas que estoy tratando de hacer:

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 = <double * > malloc((n ** 2) * sizeof(double))

            # a 2D typed memoryview onto tmp_carray
            tmp = <double[:n, :n] > 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

Parece que hay al menos dos razones por las cuales esto no se compila:

Basado en la salida decython -a, la creación de la vista de memoria escrita aquí:

cdef double[:, :] tmp = <double[:n, :n] > tmp_carray

parece involucrar llamadas API de Python, y por lo tanto no puedo liberar el GIL para permitir que el bucle externo se ejecute en paralelo.

Tenía la impresión de que las vistas de memoria escrita no eran objetos de Python y, por lo tanto, un proceso secundario debería poder crearlas sin adquirir primero el GIL. ¿Es este el caso?

2. Incluso si reemplazoprange(m, nogil=True) con un normalrange(m), A Cython todavía no parece gustarle la presencia de uncdef dentro del bucle interno:

    Error compiling Cython file:
    ------------------------------------------------------------
    ...
                # allocate a temporary array to hold the result of
                # expensive_function_1
                tmp_carray = <double*> malloc((n ** 2) * sizeof(double))

                # a 2D typed memoryview onto tmp_carray
                cdef double[:, :] tmp = <double[:n, :n]> tmp_carray
                    ^
    ------------------------------------------------------------

    parallel_allocate.pyx:26:17: cdef statement not allowed here

Actualizar

Resulta que el segundo problema se resolvió fácilmente moviendo

 cdef double[:, :] tmp

fuera de lafor bucle, y simplemente asignando

 tmp = <double[:n, :n] > tmp_carray

dentro del bucle. Sin embargo, todavía no entiendo completamente por qué esto es necesario.

Ahora si trato de usarprange Llegué al siguiente error de compilación:

Error compiling Cython file:
------------------------------------------------------------
...
            # allocate a temporary array to hold the result of
            # expensive_function_1
            tmp_carray = <double*> malloc((n ** 2) * sizeof(double))

            # a 2D typed memoryview onto tmp_carray
            tmp = <double[:n, :n]> tmp_carray
               ^
------------------------------------------------------------

parallel_allocate.pyx:28:16: Memoryview slices can only be shared in parallel sections