numpy · seberg · Jan 9, 2025 · Jan 7, 2025 · Jan 8, 2025 · Jan 8, 2025
diff --git a/numpy/_core/src/multiarray/multiarraymodule.c b/numpy/_core/src/multiarray/multiarraymodule.c
@@ -5033,6 +5033,24 @@ PyMODINIT_FUNC PyInit__multiarray_umath(void) {
         goto err;
     }
 
+    /*
+     * Initialize the default PyDataMem_Handler capsule singleton.
+     */
+    PyDataMem_DefaultHandler = PyCapsule_New(
+            &default_handler, MEM_HANDLER_CAPSULE_NAME, NULL);
+    if (PyDataMem_DefaultHandler == NULL) {
+        goto err;
+    }
+
+    /*
+     * Initialize the context-local current handler
+     * with the default PyDataMem_Handler capsule.
+     */
+    current_handler = PyContextVar_New("current_allocator", PyDataMem_DefaultHandler);
+    if (current_handler == NULL) {
+        goto err;
+    }
+
     if (initumath(m) != 0) {
         goto err;
     }
@@ -5067,7 +5085,7 @@ PyMODINIT_FUNC PyInit__multiarray_umath(void) {
      * init_string_dtype() but that needs to happen after
      * the legacy dtypemeta classes are available.
      */
-    
+
     if (npy_cache_import_runtime(
             "numpy.dtypes", "_add_dtype_helper",
             &npy_runtime_imports._add_dtype_helper) == -1) {
@@ -5081,23 +5099,6 @@ PyMODINIT_FUNC PyInit__multiarray_umath(void) {
     }
     PyDict_SetItemString(d, "StringDType", (PyObject *)&PyArray_StringDType);
 
-    /*
-     * Initialize the default PyDataMem_Handler capsule singleton.
-     */
-    PyDataMem_DefaultHandler = PyCapsule_New(
-            &default_handler, MEM_HANDLER_CAPSULE_NAME, NULL);
-    if (PyDataMem_DefaultHandler == NULL) {
-        goto err;
-    }
-    /*
-     * Initialize the context-local current handler
-     * with the default PyDataMem_Handler capsule.
-     */
-    current_handler = PyContextVar_New("current_allocator", PyDataMem_DefaultHandler);
-    if (current_handler == NULL) {
-        goto err;
-    }
-
     // initialize static reference to a zero-like array
     npy_static_pydata.zero_pyint_like_arr = PyArray_ZEROS(
             0, NULL, NPY_DEFAULT_INT, NPY_FALSE);

diff --git a/numpy/_core/src/umath/legacy_array_method.c b/numpy/_core/src/umath/legacy_array_method.c
@@ -311,7 +311,7 @@ get_initial_from_ufunc(
         }
     }
     else if (context->descriptors[0]->type_num == NPY_OBJECT
-            && !reduction_is_empty) {
+             && !reduction_is_empty) {
         /* Allows `sum([object()])` to work, but use 0 when empty. */
         Py_DECREF(identity_obj);
         return 0;
@@ -323,13 +323,6 @@ get_initial_from_ufunc(
         return -1;
     }
 
-    if (PyTypeNum_ISNUMBER(context->descriptors[0]->type_num)) {
-        /* For numbers we can cache to avoid going via Python ints */
-        memcpy(context->method->legacy_initial, initial,
-               context->descriptors[0]->elsize);
-        context->method->get_reduction_initial = &copy_cached_initial;
-    }
-
     /* Reduction can use the initial value */
     return 1;
 }
@@ -427,11 +420,47 @@ PyArray_NewLegacyWrappingArrayMethod(PyUFuncObject *ufunc,
     };
 
     PyBoundArrayMethodObject *bound_res = PyArrayMethod_FromSpec_int(&spec, 1);
+
     if (bound_res == NULL) {
         return NULL;
     }
     PyArrayMethodObject *res = bound_res->method;
+
+    // set cached initial value for numeric reductions to avoid creating
+    // a python int in every reduction
+    if (PyTypeNum_ISNUMBER(bound_res->dtypes[0]->type_num) &&
+        ufunc->nin == 2 && ufunc->nout == 1) {
+
+        PyArray_Descr *descrs[3];
+
+        for (int i = 0; i < 3; i++) {
+            // only dealing with numeric legacy dtypes so this should always be
+            // valid
+            descrs[i] = bound_res->dtypes[i]->singleton;
+        }
+
+        PyArrayMethod_Context context = {
+                (PyObject *)ufunc,
+                bound_res->method,
+                descrs,
+        };
+
+        int ret = get_initial_from_ufunc(&context, 0, context.method->legacy_initial);
+
+        if (ret < 0) {
+            Py_DECREF(bound_res);
+            return NULL;
+        }
+
+        // only use the cached initial value if it's valid
+        if (ret > 0) {
+            context.method->get_reduction_initial = &copy_cached_initial;
+        }
+    }
+
+
     Py_INCREF(res);
     Py_DECREF(bound_res);
+
     return res;
 }
diff --git a/numpy/_core/tests/test_multithreading.py b/numpy/_core/tests/test_multithreading.py
@@ -120,3 +120,17 @@ def legacy_125():
 
     task1.start()
     task2.start()
+
+def test_parallel_reduction():
+    # gh-28041
+    NUM_THREADS = 50
+
+    b = threading.Barrier(NUM_THREADS)
+
+    x = np.arange(1000)
+
+    def closure():
+        b.wait()
+        np.sum(x)
+
+    run_threaded(closure, NUM_THREADS, max_workers=NUM_THREADS)
diff --git a/numpy/testing/_private/utils.py b/numpy/testing/_private/utils.py
@@ -2685,9 +2685,9 @@ def _get_glibc_version():
 _glibc_older_than = lambda x: (_glibcver != '0.0' and _glibcver < x)
 
 
-def run_threaded(func, iters, pass_count=False):
+def run_threaded(func, iters, pass_count=False, max_workers=8):
     """Runs a function many times in parallel"""
-    with concurrent.futures.ThreadPoolExecutor(max_workers=8) as tpe:
+    with concurrent.futures.ThreadPoolExecutor(max_workers=max_workers) as tpe:
         if pass_count:
             futures = [tpe.submit(func, i) for i in range(iters)]
         else: