Repository: cython/cython
Pull request: Fix ufunc GIL handling
It looks like Numpy releases the GIL for us (unless it's PyObject arguments), so we should assume that we don't have the GIL, and write code to regain it if needed, rather than assuming we need to release the GIL.

Fixes #5328

Problem statement:
Fatal Python error: PyEval_SaveThread when creating ufunc with nogil [BUG] 
### Describe the bug

I'm generating ufuncs for my fortran library by the new version 3.0.0b1 of cython, using the @cython.ufunc function decorator.
I declare the function as nogil. and call it in python script. It works fine when the input array length is less than 500. When the input array length is greater than 500, it stops and throws a thread error.
```
Fatal Python error: PyEval_SaveThread: the function must be called with the GIL held, but the GIL is released (the current Python thread state is NULL)
Python runtime state: initialized
```
I'm using the mingw-w64 compiler provided by msys2 on windows, but that doesn't seem to be the cause of the error. I transferred the program to an ubuntu machine and found the same error.

### Code to reproduce the behaviour:

```cython
# test_ufunc.pyx
cimport cython

@cython.ufunc
cdef double add_one(double x) nogil:
    return x+1
```

```python
import test_ufunc
import numpy as np

x=np.linspace(-1,1,501)
x_add=test_ufunc.add_one(x)
print(x_add)
```

### Expected behaviour

No errors, result returns successfully.

### OS

Ubuntu and Windows

### Python version

3.10.6 (Ubuntu), 3.11.2 (Windows)

### Cython version

3.0.0b1

### Additional context

_No response_

Hints from issue discussion:
So it looks like I misunderstood ufuncs and the GIL. I'd assumed that it was up to the user to release the GIL in a ufunc (on the basis that they'd documented macros to do so https://numpy.org/doc/stable/reference/c-api/ufunc.html#macros). But it looks like Numpy releases it automatically. It only becomes an issue for large arrays though.

It also looks like I wasn't the only one to find this confusing https://github.com/numpy/numpy/issues/2067.

Should be fairly easy to fix I think.

Thanks for the report

Commit message: Fix ufunc GIL handling

It looks like Numpy releases the GIL for us (unless it's PyObject
arguments), so we should assume that we don't have the GIL, and
write code to regain it if needed, rather than assuming we need
to release the GIL.

Fixes #5328

Patch to review:
diff --git a/Cython/Compiler/UFuncs.py b/Cython/Compiler/UFuncs.py
index ce6de62fe23..2f0797bbf09 100644
--- a/Cython/Compiler/UFuncs.py
+++ b/Cython/Compiler/UFuncs.py
@@ -138,6 +138,9 @@ def generate_cy_utility_code(self):
 
         ufunc_cname = self.global_scope.next_id(self.node.entry.name + "_ufunc_def")
 
+        outer_is_nogil = not (any(t.is_pyobject for t in arg_types) or
+            any(t.is_pyobject for t in out_types))
+
         context = dict(
             func_cname=ufunc_cname,
             in_types=arg_types,
@@ -145,6 +148,7 @@ def generate_cy_utility_code(self):
             inline_func_call=self.node.entry.cname,
             inline_func_declaration=inline_func_decl,
             nogil=self.node.entry.type.nogil,
+            outer_is_nogil=outer_is_nogil,
         )
 
         code = CythonUtilityCode.load(
diff --git a/Cython/Utility/UFuncs.pyx b/Cython/Utility/UFuncs.pyx
index 6a9d3cdc652..fb2519b473d 100644
--- a/Cython/Utility/UFuncs.pyx
+++ b/Cython/Utility/UFuncs.pyx
@@ -8,7 +8,7 @@ cdef extern from *:
 
 # variable names have to come from tempita to avoid duplication
 @cname("{{func_cname}}")
-cdef void {{func_cname}}(char **args, const npy_intp *dimensions, const npy_intp* steps, void* data) except *:
+cdef void {{func_cname}}(char **args, const npy_intp *dimensions, const npy_intp* steps, void* data) except * {{"nogil" if outer_is_nogil else ""}}:
     cdef npy_intp i
     cdef npy_intp n = dimensions[0]
     {{for idx, tp in enumerate(in_types)}}
@@ -23,7 +23,7 @@ cdef void {{func_cname}}(char **args, const npy_intp *dimensions, const npy_intp
     cdef npy_intp step_{{idx}} = steps[{{idx}}]
     {{endfor}}
 
-    with nogil({{nogil}}):
+    {{"with gil" if (not nogil and outer_is_nogil) else "if True"}}:
         for i in range(n):
             {{for idx, tp in enumerate(in_types)}}
             {{if tp.is_pyobject}}
diff --git a/tests/run/ufunc.pyx b/tests/run/ufunc.pyx
index 8c9c1eefc94..7d1d3064afc 100644
--- a/tests/run/ufunc.pyx
+++ b/tests/run/ufunc.pyx
@@ -10,6 +10,17 @@ int_arr_1d = np.arange(20, dtype=int)[::4]
 int_arr_2d = np.arange(500, dtype=int).reshape((50, -1))[5:8, 6:8]
 double_arr_1d = int_arr_1d.astype(np.double)
 double_arr_2d = int_arr_2d.astype(np.double)
+# Numpy has a cutoff at about 500 where it releases the GIL, so test some large arrays
+large_int_arr_1d = np.arange(1500, dtype=int)
+large_int_arr_2d = np.arange(1500*600, dtype=int).reshape((1500, -1))
+large_double_arr_1d = large_int_arr_1d.astype(np.double)
+large_double_arr_2d = large_int_arr_2d.astype(np.double)
+
+def dont_print(x):
+    """
+    For the large arrays, it's just a "don't crash" test
+    """
+    pass
 
 # it's fairly hard to test that nogil results in the GIL actually
 # being released unfortunately
@@ -29,6 +40,8 @@ def test_tripple_it():
     array([[168., 171.],
            [198., 201.],
            [228., 231.]])
+    >>> dont_print(tripple_it(large_int_arr_1d))
+    >>> dont_print(tripple_it(large_int_arr_2d))
     """
 
 @cython.ufunc
@@ -41,6 +54,8 @@ def test_to_the_power():
     True
     >>> np.allclose(to_the_power(1., double_arr_2d), np.ones_like(double_arr_2d))
     True
+    >>> dont_print(to_the_power(large_double_arr_1d, -large_double_arr_1d))
+    >>> dont_print(to_the_power(large_double_arr_2d, -large_double_arr_2d))
     """
 
 @cython.ufunc
@@ -56,6 +71,7 @@ def test_py_return_value():
     >>> py_return_value(double_arr_1d).dtype
     dtype('O')
     >>> py_return_value(-1.)  # returns None
+    >>> dont_print(py_return_value(large_double_arr_1d))
     """
 
 @cython.ufunc
@@ -66,6 +82,7 @@ def test_py_arg():
     """
     >>> py_arg(np.array([1, "2.0", 3.0], dtype=object))
     array([1., 2., 3.])
+    >>> dont_print(py_arg(np.array([1]*1200, dtype=object)))
     """
 
 @cython.ufunc
@@ -156,3 +173,31 @@ def test_nested_function():
     >>> nested_function(-1.)
     -2.0
     """
+
+@cython.ufunc
+cdef double can_throw(double x):
+    if x<0:
+        raise RuntimeError
+    return x
+
+def test_can_throw():
+    """
+    >>> arr = double_arr_1d.copy()
+    >>> arr[1] = -1.
+    >>> can_throw(arr)
+    Traceback (most recent call last):
+    ...
+    RuntimeError
+    >>> large_arr = large_double_arr_1d.copy()
+    >>> large_arr[-4] = -2.
+    >>> can_throw(large_arr)
+    Traceback (most recent call last):
+    ...
+    RuntimeError
+    >>> large_arr2d = large_double_arr_2d.copy()
+    >>> large_arr2d[100, 200] = -1.
+    >>> can_throw(large_arr2d)
+    Traceback (most recent call last):
+    ...
+    RuntimeError
+    """


Code review comments:
diff_hunk: @@ -29,6 +40,8 @@ def test_tripple_it():
     array([[168., 171.],
            [198., 201.],
            [228., 231.]])
+    >>> dont_print(tripple_it(large_int_arr_1d))
+    >>> dont_print(tripple_it(large_int_arr_2d))
path: tests/run/ufunc.pyx
text: @user1:
It's common to use `_ = something()` in Python for (return) values that are explicitly ignored. Not sure if that's better, but it seems less verbose.

@author:
Done. Also I'd spelt `triple` wrong (in the initial commit)