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 equal the previous loop's for the layouts this branch accepts.
• The flat loops live in an out-of-line NPY_NOINLINE helper, contig_<kind>_<TYPE>, so the strided BINARY_LOOP in the caller compiles exactly as on main. Inlining them changed its register allocation and slowed strided and scalar-broadcast operands by a few percent.
• out is in1, out is in2, and out is in1 is in2 (z += z) each get their own loop. With a single three-pointer loop, icx's runtime alias check sees out == in as an overlap and runs its scalar fallback, which made in-place calls slower than main; gcc and clang vectorize the three-pointer loop either way.
• 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 -Wpass-failed warnings for CDOUBLE. The icx warning count is unchanged from main (99, of which 91 are -Wpass-failed).

Codegen (icx 2026.1.1, conda-forge CFLAGS -march=nocona -mtune=haswell -O2), from the disassembly of libmkl_umath_loops.so:

• complex64: all four loops are vectorized at 128 bits (addps/subps), with a runtime dependence check on the two-pointer loops.
• complex128: not loop-vectorized; icx unrolls and SLP-vectorizes the loops instead (addpd/subpd). The complex128 gains below come from that.

Benchmarks

Time relative to main, so below 1 is faster. Ranges span the sizes listed and both complex dtypes, for add and subtract, calling mkl_umath._ufuncs directly.

Setup: Core i7-1185G7 (Tiger Lake) under WSL2, conda-forge icx 2026.1.1, NumPy 2.5.3, MKL 2026.1.0 (gnu_thread), conda-forge default CFLAGS, MKL_NUM_THREADS=1, pinned to one core with taskset. Three builds (main, the previous revision 2bbe407, this revision) were run interleaved over 8 rounds; each number is the per-build minimum over rounds (min of 9 repeats each) divided by the main minimum. A duplicate main build run the same way gives the noise floor: 0.85–1.14 across all cases, 0.94–1.08 for 90% of them. WSL2 does not expose the CPU governor, so turbo stayed on.

Path	2bbe407	This revision
out=c (disjoint), n = 512–100000	0.46–0.93	0.45–0.89
out=a / out=b (in-place), n = 512–100000	0.45–0.86	0.43–0.89
f(a, a, out=c), n = 512–100000	0.44–0.88	0.46–0.91
f(a, a, out=a) (z += z), n = 512–100000	0.95–1.89	0.37–0.79
All modes, n = 8 and 64	0.90–1.10	0.88–1.08
Strided a[::2] operands (unchanged BINARY_LOOP), n = 8–100000	0.83–1.09	0.90–1.05
Scalar-broadcast operand (unchanged BINARY_LOOP), n = 8–100000	0.90–1.16	0.89–1.10
oneMKL VM path, n = 100001	0.86–1.02	0.96–1.10
accumulate, n = 1000	0.97–1.04	0.94–1.04

By dtype for this revision at n = 512–100000: complex64 out=c 0.45–0.89, in-place 0.43–0.89, z += z 0.52–0.79; complex128 out=c 0.60–0.89, in-place 0.55–0.84, z += z 0.37–0.73. The z += z slowdown of 2bbe407 was complex64 only (1.20–1.89); complex128 was 0.95–1.15. The strided, scalar-broadcast, oneMKL and accumulate rows for this revision are inside the noise floor. The few-percent regression on strided and scalar-broadcast operands that @antonwolfy measured for 2bbe407 is below what this host can resolve.

Tests

test_complex_contig covers out-of-place, out=a and out=b at n = 1, 7, 100, 1001 and 10**5 (the largest size this branch handles) with exact value equality against NumPy (assert_array_equal). test_complex_contig_same_input covers in1 is in2, both out-of-place and z += z. test_complex_contig_signed_zero checks every (±0, ±0) pair by sign bit, since assert_array_equal treats -0.0 and 0.0 as equal. test_accumulate_complex covers the partial-overlap path; it fails on a GCC build of an earlier revision of this PR (acced6b) and passes with this one.

The full suite passes on builds with icx 2026.1.1 and gcc 16.2 (conda-forge).

[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]

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]

Re-measured at 1346ce2 (e5bf8d2 plus the merge of main) against main 8ca254c on a second Lunar Lake host. The speedups hold, nothing outside the new branch got slower, and the only cost left is on tiny complex64 arrays.

Setup. Core Ultra 5 236V, Ubuntu 25.10. conda-forge icx 2026.1.1, NumPy 2.5.4, MKL 2026.1.0 (gnu_thread), conda-forge default CFLAGS (-march=nocona -mtune=haswell -O2 …). Pinned to P-core 2 with the performance governor and turbo off. Three builds were interleaved: main, a second main build (byte-identical, as an A/A control) and this PR. Inner loop: 41 rounds, with each build's libmkl_umath_loops.so loaded with dlopen into one process and called alternately on the same buffers. Through NumPy (mkl_umath._ufuncs): 31 rounds, one process per build, equal-length site paths. Results are medians over rounds with bootstrap 95% CIs; the A/A control stayed within ±0.5% for 95% of cases. A repeat pass at default turbo gave the same verdicts. The full suite passes (440 tests).

Gains. Time relative to main, so below 1 is faster; n = 512–100000.

Path	NumPy, complex64	NumPy, complex128	Inner loop, both dtypes
out=c	0.59–0.78	0.71–0.96	0.41–0.99
out=a / out=b	0.57–0.79	0.68–0.87	0.41–0.92
f(a, a, out=c)	0.61–0.78	0.70–0.92	0.41–0.94
f(a, a, out=a) (z += z)	0.35–0.76	0.46–0.71	0.32–0.63

At n = 64, complex128 is 0.96–0.99 through NumPy and 0.45–0.59 in the inner loop.

Unchanged paths.

Path	NumPy	Inner loop
Strided a[::2] operands, n = 8–100000	0.99–1.01	0.94–1.04
Scalar-broadcast operand, n = 8–100000	0.99–1.01	0.94–1.05
a - scalar (new output), n = 8–100000	0.99–1.01
oneMKL path, n = 100001 and 1e6, 1 thread	1.00–1.01	1.00–1.01
oneMKL path, 4 threads on P-cores 0–3	0.98–1.01
accumulate / reduce, n = 1000	1.00–1.01	1.00
Other loops: complex multiply, real add / multiply, sin, exp, absolute	1.00	1.00
npbench mandelbrot2 kernel with patch_numpy_umath(), M / L	0.976 / 1.00

The inner-loop spread on strided and scalar operands (at most +1.4% at n ≥ 64, +0.5 ns at n = 8) is code placement: the decoy build below moves unrelated cases by up to ±2% as well.

What's left: complex64 at n = 8–64 through NumPy. These calls are 1.00–1.04 against main, up to +13 ns on a ~330 ns call; the worst is subtract with out=b at n = 8. It shows even at n = 64, where the inner loop alone is 11 ns faster (33.6 → 22.6 ns), so it's a per-call cost outside the loop body. In the 2026-10-14 run it was inside the ±3% noise band. To separate code placement from the new code path, I built a decoy: this PR with the helper call behind a condition that's never true (n < 0), so it has the same layout and caller but always runs BINARY_LOOP. On the same cases the decoy was 1.00–1.01 where this PR was 1.01–1.04, so most of the cost comes from running the helper. That host has no perf counters, so I can't say which front-end effect it is.

I don't think this needs a change. It's at most 4% on calls that are almost all NumPy overhead, the gain starts between n = 64 and n = 512, and a size threshold would still keep ~1% from the layout change while giving up the complex128 gain at n = 64. It may be worth a row in the description's table.

[antonwolfy]

Lets wait for #292 merged first to have the ASV benchmark results