perf: add contiguous fast path to complex add/subtract fallback loop

[intel-python-devops]

Note

This pull request is AI-generated. Please review thoroughly.

Adds a contiguous fast path to the complex add/subtract loops (CFLOAT, CDOUBLE) for calls that don't take the oneMKL VM path, i.e. contiguous arrays with n <= VML_ASM_THRESHOLD (100000). Before this change those calls fell through to the strided BINARY_LOOP.

• The branch is gated on can_vectorize: contiguous, and each input either disjoint from the output or exactly the output. Partially overlapping operands, such as accumulate/cumsum (out == in1 + 1), still take BINARY_LOOP.
• Real and imaginary parts are processed as one flat run of 2*n scalars. Since +/- is componentwise and each element is computed independently, results are identical to the previous loop for the layouts this branch accepts.
• out is in1 and out is in2 get their own two-pointer loops. With a single three-pointer loop, the compiler's runtime alias check sees out == in as an overlap and runs its scalar fallback, which made in-place calls slower than main (icx, complex64: 0.54–0.66x at n = 1k–100k).
• No vectorization pragma. NPY_PRAGMA_VECTOR is #pragma GCC ivdep under GCC, which is unsafe for partial overlap. Under icx it did not change codegen and added two -Wpass-failed warnings for CDOUBLE. The icx warning count is now the same as main (91).

Codegen (icx 2026.1.1, -qopt-report=3, conda-forge CFLAGS -march=nocona -mtune=haswell -O2):

• complex64: all three loops are vectorized (VL 4), with a runtime dependence check.
• complex128: not loop-vectorized ("vectorization possible but seems inefficient" at 128-bit width); icx unrolls and SLP-vectorizes it instead. The complex128 gains below come from that.

Tests: test_complex_contig covers out-of-place, out=a and out=b at n = 1, 7, 100, 1001 with a bitwise comparison to NumPy. test_accumulate_complex covers the partial-overlap path; it fails on a GCC build of the previous revision of this PR (acced6b) and passes with this one.

[jharlow-intel]

unsure if auto-generated PRs like this should be putting these changes into the CHANGELOG or not. Thoughts?

[ndgrigorian]

unsure if auto-generated PRs like this should be putting these changes into the CHANGELOG or not. Thoughts?

Tough call, I would say ideally yes but I would want to make sure we have some kind of anti-AI-slop-text skill or guard before that

[antonwolfy]

We need a real benchmark numbers in the PR before it merges, including both out-of-place with no overlap and in-place paths

[jharlow-intel]

Added to the description. In-place numbers showed a regression on icx (0.54–0.66x), so in-place calls now get their own loops (1.2–2.1x).

[antonwolfy]

Requesting changes. The new loops are correct, but measurements on a Lunar Lake host show two slowdowns against main, and the benchmark table requested earlier is no longer in the description. The inline comment on mkl_umath_loops.c.src has a patch that removes both slowdowns; it was measured below as V3b.

Correctness: no issues. This PR passes the full suite (412 tests), bitwise comparisons against NumPy covering every aliasing combination, n = 0 to 200003, ±0, inf and NaN, and an accumulate/cumsum check. That holds on icx, gcc, gcc -mavx2 and clang 23. The original GCC accumulate bug is fixed, and test_accumulate_complex fails on a GCC build of acced6b, as stated.

Setup. Core Ultra 5 236V, Lunar Lake, Ubuntu 25.10. conda-forge icx 2026.1.1, NumPy 2.3.5, MKL 2026.1.0 (gnu_thread). Every build uses the same command with conda-build's default CFLAGS (-march=nocona -mtune=haswell -O2 …). Runs are pinned to P-core 2 with the performance governor and turbo off, and builds are interleaved: 21 reps for the inner loop (each build's libmkl_umath_loops.so loaded with dlopen) and 11 reps through NumPy (equal-length site paths). Results are bootstrap 95% CIs, judged against a main-vs-main noise band of ±3%. A confirmation pass at the default turbo settings gave the same verdicts.

Results. Time relative to main, so below 1 is faster. Ranges span sizes and both complex dtypes, and bold marks slower than main beyond the noise band.

• V2 = this PR at 2bbe407.
• V3a = V2 + Fix A.
• V3b = V2 + Fix A + Fix B, the inline patch.

Path	V2	V3a	V3b
out=c, f(a,a,out=c), NumPy, n ≥ 512	0.60–0.99	0.60–0.99	0.60–0.99
out=a / out=b, NumPy, n ≥ 512	0.58–0.94	0.58–0.88	0.58–0.95
out=c, inner loop, n ≥ 512	0.40–1.01	0.40–1.01	0.41–1.01
All modes except z += z, NumPy, n ≤ 64	0.96–1.02	0.96–1.02	0.96–1.03
f(a,a,out=a) (z += z), NumPy, n ≥ 512	1.23–1.86	0.41–0.85	0.41–0.86
f(a,a,out=a), NumPy, n = 64	1.12–1.15	0.97–1.02	0.97–1.03
f(a,a,out=a), inner loop, n ≥ 512	0.99–1.68	0.32–0.50	0.32–0.50
Strided / scalar-broadcast operands (unchanged BINARY_LOOP), inner loop, n = 64–1000	0.95–1.13	0.95–1.13	0.96–1.02
Same, n = 8	0.94–1.19	0.96–1.19	0.96–1.06 (+0.3–0.7 ns)
NumPy complex64 a - scalar, n = 1000	1.035 [1.031, 1.038]	1.077 [1.069, 1.078]	1.015 [1.009, 1.026]
oneMKL path, n = 100001 and 1e6, 1 and 4 threads	0.98–1.02	0.98–1.02	0.98–1.02
accumulate / cumsum, n = 1000	1.00	1.00	1.00
npbench mandelbrot2 M / L	0.99 / 1.00	0.99 / 1.00	0.99 / 1.00

An earlier run on V2 measured further legacy paths, all within noise: complex multiply, add/subtract reduce, real float32/float64 add/subtract/multiply/divide, sin/exp/log/sqrt, the full ASV suite (174 benchmarks) and npbench.

Causes. Both are visible in the disassembly.

1. z += z. All three pointers are equal, so the out == in1 loop reads through two pointers that alias. icx's runtime overlap check fails, and it runs its scalar fallback one element at a time. Fix A gives this case a single-pointer loop, which vectorizes (addps/addpd) with no overlap check.
2. Strided / scalar-broadcast operands. The BINARY_LOOP in the same functions wasn't edited, but adding the flat loops changed its register allocation. Fix B moves the flat loops into a NPY_NOINLINE helper, which restores main's code for that loop. Instruction counts for the loop body / remainder loop:

Function	main	V2	V3a	V3b
CFLOAT_add	27 / 15	29 / 15	29 / 15	27 / 15
CFLOAT_subtract	29 / 15	35 / 24 (+imul)	35 / 24	27 / 15
CDOUBLE_add	33 / 15	35 / 15	35 / 15	33 / 15
CDOUBLE_subtract	35 / 15	41 / 24 (+imul)	41 / 24	33 / 15

What V3b still costs.

• In the bare inner loop at n ≤ 8, +0.3–2.5 ns per call for the extra function call. This doesn't show through NumPy.
• complex64 subtract out=b at n = 512–4096 is up to 18% behind V2's speedup through NumPy, though still 0.60–0.95× main. The vector loops are instruction-identical to V2's, so this looks like a code-placement effect.

The dropped benchmark table, checked. It used speedup (higher is better):

Claim	Measured
In-place out=a / out=b 1.2–2.1x at n = 1024 / 16384 / 100000	Confirmed: 1.13–1.60x through NumPy, 1.17–2.41x inner loop. complex128 at n = 100000 is 1.13–1.18x.
complex64 add out-of-place 1.48x at n = 1024	Confirmed: 1.49x
complex64 add out-of-place 1.94x at n = 16384	Not reproduced: 1.54–1.61x on this host
~1.0x at n = 64	Confirmed
mandelbrot2 1.01x	Confirmed

Requested changes

1. Apply the inline patch (V3b), or an equivalent that removes both slowdowns. Add z += z and strided/scalar-broadcast rows to the benchmark table.
2. Restore the Benchmarks section of the description. It was added in the 2026-10-09 18:14 edit and removed in the 18:42 edit, so "The complex128 gains below" now refers to nothing. Also, "bitwise comparison" is really exact value equality.
3. Tests: see the inline suggestions for in1 is in2 coverage (including z += z), n = 10**5 (the largest size this branch handles), and an optional signed-zero test.

FYI, outside this PR:

• No CI job builds with GCC, so test_accumulate_complex can only catch the original bug on a local GCC build.
• The GitHub Actions test_windows jobs run zero tests. conda activate without call ends the cmd step before pytest starts; see the test_windows (3.13) log on this PR.

[antonwolfy]

Patch for both slowdowns (V3b in the review summary). GitHub can't make this a one-click suggestion, because the helper has to be defined above the function, outside this PR's diff. It applies on top of 2bbe407 with git apply:

diff --git a/mkl_umath/src/mkl_umath_loops.c.src b/mkl_umath/src/mkl_umath_loops.c.src
index 16aea16f..a1ae2da0 100644
--- a/mkl_umath/src/mkl_umath_loops.c.src
+++ b/mkl_umath/src/mkl_umath_loops.c.src
@@ -1282,6 +1282,63 @@ pairwise_sum_@TYPE@(@ftype@ *rr, @ftype@ * ri, char * a, npy_intp n, npy_intp st
  * #PW = 1, 0#
  * #VML = Add, Sub#
  */
+/*
+ * Kept out of line so that adding these loops does not change the
+ * register allocation of the strided BINARY_LOOP in the caller.
+ */
+NPY_NOINLINE void
+contig_@kind@_@TYPE@(char **args, npy_intp n)
+{
+    /*
+     * Contiguous and either disjoint or exactly in-place, but too
+     * small for the VML call in the caller. +/- acts on the real and
+     * imaginary parts independently, so treat them as one flat
+     * run of 2*n @ftype@ elements. The in-place cases get their
+     * own loops so the compiler's runtime alias check compares
+     * only the two distinct buffers; with three pointers it sees
+     * out == in as an overlap and runs the scalar fallback.
+     * Partially overlapping operands (e.g. accumulate, where
+     * out == in1 + 1) take the caller's strided loop.
+     */
+    const npy_intp n2 = n * 2;
+    npy_intp i;
+
+    if (args[2] == args[0]) {
+        @ftype@ *iop1 = (@ftype@ *)args[0];
+
+        if (args[1] == args[0]) {
+            /* z += z: a single pointer leaves no alias check to fail */
+            for (i = 0; i < n2; i++) {
+                iop1[i] = iop1[i] @OP@ iop1[i];
+            }
+        }
+        else {
+            const @ftype@ *ip2 = (const @ftype@ *)args[1];
+
+            for (i = 0; i < n2; i++) {
+                iop1[i] = iop1[i] @OP@ ip2[i];
+            }
+        }
+    }
+    else if (args[2] == args[1]) {
+        const @ftype@ *ip1 = (const @ftype@ *)args[0];
+        @ftype@ *iop2 = (@ftype@ *)args[1];
+
+        for (i = 0; i < n2; i++) {
+            iop2[i] = ip1[i] @OP@ iop2[i];
+        }
+    }
+    else {
+        const @ftype@ *ip1 = (const @ftype@ *)args[0];
+        const @ftype@ *ip2 = (const @ftype@ *)args[1];
+        @ftype@ *op1 = (@ftype@ *)args[2];
+
+        for (i = 0; i < n2; i++) {
+            op1[i] = ip1[i] @OP@ ip2[i];
+        }
+    }
+}
+
 void
 mkl_umath_@TYPE@_@kind@(char **args, const npy_intp *dimensions, const npy_intp *steps, void *NPY_UNUSED(func))
 {
@@ -1307,45 +1364,7 @@ mkl_umath_@TYPE@_@kind@(char **args, const npy_intp *dimensions, const npy_intp
             return;
         }
         else if (can_vectorize) {
-            /*
-             * Contiguous and either disjoint or exactly in-place, but too
-             * small for the VML call above. +/- acts on the real and
-             * imaginary parts independently, so treat them as one flat
-             * run of 2*n @ftype@ elements. The in-place cases get their
-             * own loops so the compiler's runtime alias check compares
-             * only the two distinct buffers; with three pointers it sees
-             * out == in as an overlap and runs the scalar fallback.
-             * Partially overlapping operands (e.g. accumulate, where
-             * out == in1 + 1) take the strided loop below.
-             */
-            const npy_intp n2 = dimensions[0] * 2;
-            npy_intp i;
-
-            if (args[2] == args[0]) {
-                @ftype@ *iop1 = (@ftype@ *)args[0];
-                const @ftype@ *ip2 = (const @ftype@ *)args[1];
-
-                for (i = 0; i < n2; i++) {
-                    iop1[i] = iop1[i] @OP@ ip2[i];
-                }
-            }
-            else if (args[2] == args[1]) {
-                const @ftype@ *ip1 = (const @ftype@ *)args[0];
-                @ftype@ *iop2 = (@ftype@ *)args[1];
-
-                for (i = 0; i < n2; i++) {
-                    iop2[i] = ip1[i] @OP@ iop2[i];
-                }
-            }
-            else {
-                const @ftype@ *ip1 = (const @ftype@ *)args[0];
-                const @ftype@ *ip2 = (const @ftype@ *)args[1];
-                @ftype@ *op1 = (@ftype@ *)args[2];
-
-                for (i = 0; i < n2; i++) {
-                    op1[i] = ip1[i] @OP@ ip2[i];
-                }
-            }
+            contig_@kind@_@TYPE@(args, dimensions[0]);
         }
         else {
             BINARY_LOOP {

• Fix A, the args[1] == args[0] branch: z += z gets a single-pointer loop. icx has no aliasing to check there, so it vectorizes instead of falling back to scalar. Subtract still computes x - x, so NaN/inf inputs give NaN; this was checked bitwise.
• Fix B, the NPY_NOINLINE helper: the strided BINARY_LOOP in the caller compiles as it does on main again (instruction counts in the summary). NPY_NOINLINE is defined in numpy/npy_common.h as static plus noinline, with an MSVC branch, in NumPy 1.26.4 and 2.x.
• Verified with icx 2026.1.1 and gcc 15.2 (-mavx2): full suite (412 passed), the suggested tests, bitwise checks, and accumulate/cumsum for n = 2..300. Only the four complex add/subtract functions change; the vector loops are instruction-identical to the current revision's.

Wording nit in the moved comment: "the compiler's runtime alias check … runs the scalar fallback" describes icx only. gcc 13.3 and clang 23 vectorize the three-pointer loop with out == in.

[antonwolfy]

n = 10**5 is the largest size this branch handles (the oneMKL call starts at n > 100000), so it's worth covering:

Suggested change

	@pytest.mark.parametrize("size", [1, 7, 100, 1001])
	@pytest.mark.parametrize("size", [1, 7, 100, 1001, 10**5])

[antonwolfy]

assert_array_equal checks value equality: it treats -0.0 == +0.0 and only compares NaN positions. FP add also rounds, so "+/- is exact" doesn't hold.

Suggested change

	# and in-place on either input; +/- is exact, so compare bitwise
	# and in-place on either input; results must equal NumPy's exactly

[antonwolfy]

Coverage for in1 is in2, both out-of-place and z += z (each takes a different loop), plus an optional signed-zero test, since assert_array_equal can't see the sign of zero. Both tests pass on main, on this revision and with the inline patch, on icx and gcc.

Suggested change

	assert_array_equal(mkl_res, np_res)
	assert_array_equal(mkl_res, np_res)
	@pytest.mark.parametrize("func", ["add", "subtract"])
	@pytest.mark.parametrize("size", [1, 7, 100, 1001, 10**5])
	@pytest.mark.parametrize("dtype", [np.complex64, np.complex128])
	def test_complex_contig_same_input(func, size, dtype):
	# in1 is in2: out-of-place takes the disjoint loop with ip1 == ip2,
	# out=a takes the out == in1 loop with ip2 == iop1
	a = _complex_array(size, dtype)
	mkl_func = getattr(mu, func)
	np_res = getattr(np, func)(a, a)
	assert_array_equal(mkl_func(a, a), np_res)
	a_inplace = a.copy()
	mkl_func(a_inplace, a_inplace, out=a_inplace)
	assert_array_equal(a_inplace, np_res)
	@pytest.mark.parametrize("func", ["add", "subtract"])
	@pytest.mark.parametrize("dtype", [np.complex64, np.complex128])
	def test_complex_contig_signed_zero(func, dtype):
	# every (+-0, +-0) pair in both the real and the imaginary part;
	# assert_array_equal treats -0.0 == 0.0, so compare the sign bits
	ftype = np.float32 if dtype == np.complex64 else np.float64
	zeros = np.array([0.0, -0.0], dtype=ftype)
	x, y = np.meshgrid(zeros, zeros)
	a = np.tile(np.repeat(x.ravel(), 2), 8).view(dtype)
	b = np.tile(np.repeat(y.ravel(), 2), 8).view(dtype)
	mkl_func = getattr(mu, func)
	np_res = getattr(np, func)(a, b)
	for out_is in ("new", "a", "b"):
	a_, b_ = a.copy(), b.copy()
	out = {"new": np.empty_like(a), "a": a_, "b": b_}[out_is]
	mkl_func(a_, b_, out=out)
	assert_array_equal(
	np.signbit(out.view(ftype)),
	np.signbit(np_res.view(ftype)),
	err_msg=f"signed zero mismatch for '{func}' with out={out_is}",
	)