Optimize Standard complex CPU contractions with faer GEMM

[exAClior]

Summary

Standard<Complex32> and Standard<Complex64> CPU contractions currently miss the faer GEMM fast paths used by Standard<f32> and Standard<f64>. They therefore fall through to generic_gemm, the scalar triple loop, even when contraction planning has already lowered the operation to an ordinary GEMM.

This issue tracks the benchmark, implementation, correctness proof, and downstream validation needed to merge optimized complex dispatch without adding a TDVP-specific workaround to OMEinsum.

The motivating downstream workload is QudeLeap/rydbergsim-rs#272, where complex local-Hamiltonian contractions dominate saturated two-site TDVP runtime.

Current diagnosis

At main commit 239e962, optimized CPU dispatch is limited to real standard algebra in four entry points in src/backend/cpu/mod.rs:

1. Cpu::gemm_internal
2. Cpu::gemm_standard_layout_internal
3. Cpu::gemm_batched_internal
4. Cpu::gemm_batched_standard_layout_internal

The contraction planner already classifies modes, handles directly usable layouts, materializes required permutations, and lowers binary contractions to GEMM. The missing piece is complex scalar dispatch at that GEMM boundary.

Work prepared

A benchmark-and-plan branch exists at:

• Branch: exAClior/omeinsum-rs:perf/complex-gemm
• Scaffold commit: d6695fc
• Detailed plan: benchmarks/complex_tdvp.md
• Criterion target: benches/complex_tdvp.rs

The benchmark calls the public Tensor::contract_binary::<Standard<Complex64>> path with deterministic inputs, MPO bond dimension D=9, and chi in {32, 64, 128}. It covers four contractions taken directly from the downstream TDVP local-Hamiltonian applications:

Case	Equivalent GEMM
one-site left environment	(D chi × chi) @ (chi × 2 chi)
one-site right environment	(2 chi × chi D) @ (chi D × chi)
two-site left environment	(D chi × chi) @ (chi × 4 chi)
two-site right environment	(4 chi × chi D) @ (chi D × chi)

make check passes on the current comparison working tree: formatting, all-target clippy, 158 unit tests, 333 integration tests, and 16 doctests passed. The OMEinsum-versus-Tenferro results are summarized below and reported in full in this comment.

Implementation plan

• Verify that complex standard contractions reach generic_gemm.
• Add deterministic TDVP-shaped public-path benchmarks.
• Document target paths, workload, risks, and acceptance criteria.
• Record scalar baselines for chi=32, then chi=64, then chi=128 under the repository's runscribe protocol.
• Add faer dispatch for contiguous Complex32 and Complex64 GEMM.
• Add complex dispatch for borrowed strided/transposed layouts.
• Add complex dispatch for materialized and direct-layout batched GEMM.
• Add concrete-value tests against generic_gemm for both complex widths.
• Cover rectangular, transposed, positive-stride, supported negative-stride, and multi-batch layouts.
• Verify that ordinary GEMM semantics do not accidentally conjugate either operand.
• Run candidate benchmarks and existing real-valued binary benchmarks.
• Pin the candidate downstream, remove the local prototype contraction workaround, and rerun saturated TDVP kernel/step measurements.
• Open a focused upstream PR with before/after results and downstream evidence.

Implementation constraints

• Keep this inside the existing CPU GEMM boundary; no public API or contraction-order change.
• Reuse faer_mat_ref, layout bounds checking, and the existing checked-TypeId dispatch pattern.
• Borrow inputs where possible and allocate the output once.
• Keep Par::Seq initially, matching existing real-valued behavior; parallel policy is a separate benchmark question.
• Do not create a permanent downstream mini-einsum implementation.
• Prefer four short, obvious dispatch arms over a trait hierarchy unless genericization materially reduces unsafe casts or duplicated layout code.

Correctness acceptance

• Complex fast paths agree with generic_gemm within dtype-appropriate absolute and relative tolerances.
• Tests use nontrivial real and imaginary values and concrete expected values.
• Contiguous, strided/transposed, and batched paths are covered for Complex32 and Complex64.
• Existing standard-real, tropical, repeated-label, scalar-output, and backward tests remain green.
• make check passes.

Performance acceptance

On one pinned host/toolchain:

• At least 5× lower median latency than the scalar baseline for all four chi=64, D=9 target contractions.
• No target contraction regresses.
• chi=128 improves in the same direction without pathological allocation or memory growth.
• Existing real-valued binary benchmarks do not regress by more than 5% beyond measurement noise.
• Downstream saturated TDVP timing confirms that the kernel improvement survives contraction planning, environment updates, and Krylov repetition.

The 5× kernel threshold is a minimum useful result, not a claim that it alone makes the full paper-reproduction grid feasible.

OMEinsum versus Tenferro result (2026-07-15)

Both implementations resolve complex GEMM to single-threaded faer 0.24.4. With identical column-major Complex64 inputs, Tenferro's prepared contraction path ties OMEinsum on the two left-environment layouts and wins on the two right-environment layouts.

The reverse-order chi=64 check measured:

Case	OMEinsum median	Tenferro median	Tenferro speedup
h1-left-environment	0.8066 ms	0.8031 ms	1.00x
h1-right-environment	1.0950 ms	0.8784 ms	1.25x
h2-left-environment	1.6763 ms	1.6670 ms	1.01x
h2-right-environment	2.2396 ms	1.7475 ms	1.28x

At chi=128, the right-environment advantages narrow to 1.16x and 1.09x. This is not evidence for a faster Tenferro arithmetic kernel; the leading explanation is different contraction lowering, including plan reuse, strided-layout handling, materialization, buffer reuse, or per-call framework overhead. Profiling is still required to isolate the mechanism.

Full protocol, all chi values, the order-bias check, and caveats are in the comparison comment. runscribe was unavailable, so the repository owner explicitly approved a documented fallback that captured raw Criterion output, samples, confidence intervals, environment metadata, and the exact benchmark patch.

[exAClior]

OMEinsum vs Tenferro complex TDVP comparison (2026-07-15)

I added Tenferro to the existing Complex64, D=9, chi ∈ {32,64,128} benchmark and ran the same four contractions through both libraries.

Protocol

• Host: Apple M4 Pro, 12 CPUs, 64 GiB; macOS 15.7.4 arm64
• Rust: 1.93.1
• OMEinsum base: c2866c9caa4b918367002c0478d47a41050e72e4
• Tenferro: pinned to c76f260b595c292aa5fbedfc87de80b06e6a0a8f
• Identical deterministic column-major Complex64 buffers and explicit output order
• OMEinsum: Tensor::contract_binary::<Standard<Complex64>>
• Tenferro: prepared ConcreteEinsumPlan::execute_typed; plan construction excluded from timing
• Both calls return an allocated output
• Both faer paths explicitly single-threaded (Par::Seq / CpuBackend::with_threads_and_kind(1, Faer))
• Criterion: 10 samples, 2 s warmup, 10 s measurement per benchmark
• All 24 library/case/size combinations passed shape and elementwise value checks before timing
• make check passed: 158 unit, 333 integration, and 16 doctests passed

Tenferro speedup = OMEinsum median / Tenferro median; values above 1 favor Tenferro.

Primary results

Case	chi	OMEinsum median (ms)	Tenferro median (ms)	Tenferro speedup
h1-left-environment	32	0.1003	0.1003	1.00x
h1-left-environment	64	0.8143	0.8080	1.01x
h1-left-environment	128	6.7072	6.6680	1.01x
h1-right-environment	32	0.1445	0.1123	1.29x
h1-right-environment	64	1.1414	0.8910	1.28x
h1-right-environment	128	7.9358	6.8351	1.16x
h2-left-environment	32	0.1997	0.1983	1.01x
h2-left-environment	64	1.7195	1.7859	0.96x
h2-left-environment	128	12.9959	12.9277	1.01x
h2-right-environment	32	0.2970	0.2167	1.37x
h2-right-environment	64	2.2811	1.7576	1.30x
h2-right-environment	128	14.9334	13.7204	1.09x

For an illustrative workload containing one call to each contraction, this is about 1.18x at chi=32, 1.14x at chi=64, and 1.06x at chi=128. Real TDVP call frequencies may differ.

Reverse-order check at chi=64

The primary benchmark registered OMEinsum first. I reran each chi=64 implementation in a separate process with Tenferro first to test thermal/order bias:

Case	OMEinsum median (ms)	Tenferro median (ms)	Tenferro speedup
h1-left-environment	0.8066	0.8031	1.00x
h1-right-environment	1.0950	0.8784	1.25x
h2-left-environment	1.6763	1.6670	1.01x
h2-right-environment	2.2396	1.7475	1.28x

The initial apparent 4% Tenferro regression on h2-left, chi=64 did not survive reversed ordering; it was noise. The right-environment advantage did survive.

Correct conclusion and likely reason

Both libraries resolve the complex GEMM itself to single-threaded faer 0.24.4. The benchmark therefore provides no evidence that Tenferro has a faster arithmetic kernel.

The left-environment paths are tied. Tenferro's prepared public contraction path consistently wins the right-environment layouts, although the gap shrinks at chi=128. Those layouts interleave free and contracted axes, unlike the copy-friendly left-environment layouts.

The leading explanation is a difference in contraction lowering: plan reuse, strided-layout handling, materialization, buffer reuse, or other per-call framework overhead. This mechanism is plausible but not yet isolated, so it should be reported as a hypothesis rather than a proven cause. Profiling allocations and the lowered layouts for the two right-environment contractions is the useful next step.

runscribe remained unavailable; the repository owner explicitly approved a documented fallback. Raw output, Criterion samples/95% intervals, environment metadata, exact benchmark patch, and the aggregation script were captured locally under benchmarks/runs/20260715T072027Z-tenferro-vs-omeinsum/ and can be pushed with the benchmark changes.

[exAClior]

Follow-up profiling isolates the Tenferro gap to operand materialization, not a different complex GEMM.

Lowered layouts

• h1-left and h2-left: neither operand is copied; faer consumes borrowed strided matrix views.
• h1-right (abcd,bde->ace): only the left operand is packed, [a,b,c,d] -> [a,c,b,d] ([0,2,1,3]).
• h2-right (abcde,cef->abdf): only the left operand is packed, [a,b,c,d,e] -> [a,b,d,c,e] ([0,1,3,2,4]).
• No case needs an output permutation.

OMEinsum currently default-initializes a fresh scratch Vec, overwrites it with the permutation, then calls .as_slice().to_vec() and copies the complete packed payload into a second fresh Vec. Tenferro performs one blocked permutation into an uninitialized thread-local pooled buffer. The traffic table below conservatively counts only the logical copies; OMEinsum's scratch zeroing may add another payload-sized write.

At chi=64, the packed payloads are 1.125 MiB (h1-right) and 2.25 MiB (h2-right). That means:

Case	OMEinsum input-pack traffic	Tenferro input-pack traffic
h1-right	at least 4.5 MiB (two read+write copies)	at least 2.25 MiB (one read+write copy)
h2-right	at least 9 MiB	at least 4.5 MiB

Steady-state allocation deltas against a prepacked-input ablation confirm this. OMEinsum requests two payloads per call (2,359,616 / 4,718,992 extra bytes including tiny metadata); warm Tenferro requests only 384 / 432 extra metadata bytes. On a cold thread Tenferro requests one payload, then reuses it on the next call.

Timed phase split (chi=64, one thread)

Case	Library	Pack	Output allocation	faer	Pipeline
h1-right	OMEinsum	137.4 us	0.6 us	884.2 us	1,028.1 us
	Tenferro	16.3 us	0.03 us	856.1 us	872.9 us
h2-right	OMEinsum	423.7 us	5.3 us	1,744.4 us	2,195.4 us
	Tenferro	31.2 us	0.05 us	1,713.3 us	1,749.3 us

The OMEinsum pack is 8.5x / 13.6x slower. Splitting that pack internally gives initialization / scalar permutation / final clone times of 5.0 / 111.3 / 21.0 us for h1-right, and 87.5 / 202.2 / 134.0 us for h2-right. The scalar loop is expensive, but the zeroed allocation and second full copy are plainly not free. Tenferro replaces the whole sequence with one blocked copy into a warm destination.

A public-path ablation prepacked the left operand outside timing. At chi=64:

Case	Original OMEinsum / Tenferro	Prepacked OMEinsum / Tenferro	Gap removed
h1-right	1,053.6 / 870.6 us	880.0 / 856.4 us	87%
h2-right	2,145.5 / 1,735.2 us	1,724.3 / 1,711.6 us	97%

Once packing is removed, the paths are within about 0-3%, like the copy-free left cases.

One threading nuance: OMEinsum uses Par::Seq; Tenferro's no-conjugation strided path uses Par::rayon(0). RAYON_NUM_THREADS=1 was therefore required in addition to constructing Tenferro's backend with one thread. The remaining faer-phase difference is small and consistent with that policy difference; this still does not show a different or intrinsically faster complex GEMM implementation.

Follow-up scope

The generic operand-materialization optimization is deliberately out of scope for #53/#54. It is tracked separately in issue #55 and stacked draft PR #56. This comment records only the diagnosis needed to avoid misattributing the right-layout gap to complex GEMM.