perf: rework the carré du champ inner loop; retire the porting plan

[wardjm]

Two independent changes: dropping PORTING_PLAN.md now that the port is done, and a measured performance pass over the two places the time actually goes.

Retiring the porting plan

The port is complete and the repo diverges from upstream from here, so the phased plan has nothing left to track. Scope lives in README.md; the places where the Julia deliberately does not match the Python live in docs/src/upstream-bugs.md, which is where the one substantive cross-reference (the fixture-correction guard) now points.

Where the time goes

Profiling a 3000–5000 point cloud found the hot spots are not the einsums:

• from_point_cloud is ~50% Arpack eigs and ~44% markov_chain — almost all of the latter is tune_kernel's exp sweep (80 ε × 60k entries, run twice).
• Every differential operator is dominated by carre_du_champ_knn, not the OMEinsum contractions. For lie_bracket it was ~93% of the build.

carre_du_champ_knn

The per-point work now runs in a point-last layout. The neighbour gather becomes a contiguous column copy instead of a strided walk over an (n, F) array, and the rank-k update accumulates into a dense H×F buffer rather than the strided cdc[p, :, :] view — which had been forcing LinearAlgebra's generic matmul fallback once per point. Only that buffer is scattered back, and the 1/2ρ scaling moves to a single vectorised pass over the contiguous axis.

Two thresholds in here fell out of measurement rather than taste, and both are commented as such: BLAS beats a hand-rolled rank-1 accumulation only above ~2^11 flops (below that the per-point mul! call overhead dominates), and the mean-centring wants opposite loop orders for wide and narrow tensors.

Threading

The cdc point loop and the ε sweep are both embarrassingly parallel, so both are threaded. The cdc split is gated on F*H*k — ungated, spawning for small tensors cost more than it saved. Neither split changes results.

Also included: skipping exp terms that underflow to exactly zero (bit-identical, worth only a few percent), and hoisting a squared-distance matrix markov_chain was building twice.

Results

n=4000, cold geometry per rep, against a stashed baseline in the same session:

	1 thread	vs base	8 threads
grad	0.0072	−30%	0.0025
d(1)	0.0131	−20%	0.0091
laplacian(1)	0.0395	−42%	0.0241
hessian	0.1834	−10%	0.1597
levi_civita	0.0282	−8%	0.0170
lie_bracket	0.0427	−35%	0.0203
build	0.2545	~flat	0.1943

Testing

Full suite green at both 1 and 8 threads — 2473 tests, parity fixtures and doctests included. That is the real gate on the muladd and loop-reordering changes, since they perturb floating-point association.

Notes for review

• Type stability is not a lever here. The numeric kernels all infer concretely, so the usual performance-tips checklist finds nothing. The abstract fields that do exist (LinearOperator._weak::Union{Nothing,AbstractMatrix}, GammaCache::Any) are per-operator caches hit once per big matrix; tightening them buys nothing measurable, so they are untouched.
• The build is now floored by Arpack, which is external and single-threaded. Replacing it (ArnoldiMethod/KrylovKit) is the next lever if build time matters more than operator time, but it carries real reproducibility risk — the current eigs call is deliberately pinned to a fixed v0 and a tight tolerance, and the comment there explains why. Deliberately out of scope.