[kevinmoran1988]

The SSV filter is one of the largest costs in a search, and nearly all of its time is in the calc_band_N kernels. This PR makes those kernels faster on x86. Scores don't change.

What was slow

Each step of a band subtracts a score from every vector in the band and folds the result into one running maximum:

sv00 = subs(sv00, *rsc);  xEv = max(xEv, sv00);
sv01 = subs(sv01, *rsc);  xEv = max(xEv, sv01);
sv02 = subs(sv02, *rsc);  xEv = max(xEv, sv02);
...

Every max has to wait for the one before it. A band can be up to 14 vectors wide, which makes a chain of 14 operations that can only run one after another, although the subtractions feeding them are independent.

What changes

Bands of three or more vectors keep two running maxima, one for the even vectors and one for the odd ones, and join them when the band ends. Each chain is half as long, and the two run side by side. The maximum of the two is the same value the single one held, so the result is identical.

The NEON kernel already does this with six accumulators. This brings the same idea to impl_avx/ssvfilter_avx.c, impl_avx/ssvfilter_sse.c and impl_sse/ssvfilter.c. NEON is not touched.

Bands of one or two vectors keep a single maximum. There is no chain to split there, and a second one only added register setup. Those kernels compile to the same instructions as before.

Kernel timing

p7_SSVFilter() alone on one core, best of 5, change in time per call. Rows are model length, columns are sequence length.

AVX2:

model	30	100	300	1000
100	−18%	−34%	−39%	−40%
200	−29%	−38%	−41%	−42%
400	−37%	−42%	−44%	−44%
1000	−42%	−37%	−36%	−35%
2000	−45%	−43%	−38%	−36%
3841	−46%	−47%	−39%	−36%

SSE:

model	30	100	300	1000
100	−30%	−39%	−41%	−42%
200	−37%	−43%	−45%	−46%
400	−44%	−44%	−42%	−42%
1000	−45%	−43%	−35%	−31%
2000	−48%	−48%	−33%	−24%
3841	−47%	−47%	−42%	−27%

Models short enough to fit in two vectors (up to 64 positions with AVX2, 32 with SSE) are within ±2%, which is the run-to-run noise. Sequences of length 1 move by −8% to +6%, a few nanoseconds either way; the default minimum ORF length of 20 means a search doesn't make such calls.

bathsearch

Against a 352 Mb genome with 16 threads:

search	build	time	CPU time
one model, 3,841 positions	AVX2	2.14 → 1.66 s	32.1 → 24.3 s
50 models	AVX2	31.5 → 32.6 s	385 → 320 s
50 models	SSE only	38.5 → 33.5 s	543 → 420 s

The 50-model AVX2 run uses 17% less CPU but takes no less time. That run is limited by the single thread that reads the target, not by the search, so the search threads just wait longer. TravisWheelerLab/easel#6 speeds up the reader; with both changes that search takes 20.9 s (24.4 s with the reader change alone).

What I checked

• 22,000 SSV calls, on sampled models of 5 to 3,841 positions and random sequences of 1 to about 80,000 residues: return status and score are identical before and after, with AVX2 and with SSE.
• The hit tables are identical in every bathsearch run above, in 50 single-model searches of the genome, and in 523 small single-thread searches.
• make check passes on an AVX2 build and on an SSE-only build.

I have no ARM machine, which is one reason NEON is left as it is. I also couldn't test on a CPU without AVX2 other than by configuring the SSE-only build.

If you'd like this arranged differently, or would rather keep a single maximum in some of the kernels, I'm happy to change it.

🤖 Generated with Claude Code

[GenevieveKrause]

The change is correct, but we couldn't reproduce the speedup.

On 80,000 SSV calls, with models of 1 to 3,983 positions and both random and model-emitted sequences, the status and score matched main exactly in AVX2 and SSE-only builds. We timed p7_SSVFilter() alone on one pinned core, best of 5, comparing main with main plus this PR, on two machines:

• Intel Xeon Platinum 8568Y+, gcc 13.3.0
• AMD EPYC 7642, gcc 13.2.0

Change in time per call; negative means the PR is faster. Rows are model length, columns are sequence length.

Intel, AVX2:

model	30	100	300	1000
100	+1%	−4%	0%	0%
200	+3%	−3%	−4%	−4%
400	+3%	+3%	+3%	0%
1000	+22%	+3%	−1%	−3%
2000	+31%	+8%	+6%	+8%
3841	+61%	+26%	+10%	+9%

Intel, SSE:

model	30	100	300	1000
100	0%	0%	0%	−1%
200	0%	−1%	0%	0%
400	+2%	+2%	+3%	+3%
1000	+13%	+6%	+6%	+6%
2000	+21%	+14%	+6%	+6%
3841	+26%	+23%	+7%	+6%

AMD, AVX2:

model	30	100	300	1000
100	+1%	0%	−1%	−1%
200	−5%	−7%	−7%	−7%
400	−5%	−12%	−15%	−15%
1000	+26%	−7%	−15%	−18%
2000	+30%	+5%	−1%	−4%
3841	+37%	+21%	+7%	+5%

AMD, SSE:

model	30	100	300	1000
100	−5%	−9%	−9%	−10%
200	0%	−10%	−17%	−17%
400	+7%	+1%	0%	−1%
1000	+20%	+8%	+3%	+3%
2000	+28%	+19%	+10%	+6%
3841	+26%	+20%	+6%	+5%

On the AMD machine, mid-size models run up to 18% faster. Neither machine comes close to the 35–47% in your
tables. Large models with short sequences are slower on both: 13–61% slower at L = 30 for models of 1,000
positions or more.

Which CPU and compiler did you use, and how did you time the kernel (build flags, how the benchmark called
p7_SSVFilter(), how many iterations)?

[kevinmoran1988]

You're right, and thank you for testing this so carefully. The gain in my tables is specific to my CPU, and the slowdown you measured for long models is real. I've reworked the kernel so that it no longer needs the extra register behind the slowdown. Details below.

Your questions

• CPU: AMD Ryzen 9 9950X (Zen 5).
• Compiler and flags: gcc 15.2.0 with the default ./configure, so -O3, plus -mavx2 for the AVX file.
• Benchmark: a model from p7_hmm_Sample(), 64 random sequences per cell drawn at the background frequencies, 2,000,000 / (L × (M/32 + 1)) + 20 calls per pass, best of 5 passes. It calls p7_SSVFilter_avx() and p7_SSVFilter_sse() directly. The core wasn't pinned for the tables in the description; pinned, I get the same numbers. I can post the harness if you'd like to compare it with yours.

Why you saw no gain

The change shortens a chain of max operations, and that only helps where the CPU is waiting on the chain. Take a band of 13 vectors. Each step does 13 saturating subtracts (_mm_subs_epi8), one per vector, and 13 byte maxes (_mm_max_epu8). Two things limit how fast a step can go:

• Units. A core has a few execution units for vector integer work. Each starts one operation per cycle, and not every unit does every operation.
• The chain. In main every max takes the result of the max before it, 13 in a row, and the first takes the last one of the step before.

What the CPUs can do:

CPU	byte max: units	byte max: cycles when chained	saturating subtract: units	saturating subtract: cycles when chained
Intel Emerald Rapids	2	1	the same 2	1
Zen 2	3	1	2	1
Zen 3	4	1	2	1
Zen 4	4	1	2	2
Zen 5	4	2	2	2

A step takes the larger of what the units need for its 26 operations and what the chain needs. In cycles:

CPU	units need	chain, main	chain, PR as posted	chain, fold (below)
Intel Emerald Rapids	13	13	7	1
Zen 2	9	13	7	1
Zen 3	7	13	7	1
Zen 4	7	13	7	2
Zen 5	7	26	14	2

• On Intel the units need 13 cycles however short the chain is, so there was nothing to gain.
• On Zen 2 a step can go from 13 cycles to 9, which fits the gain you saw on mid-size models.
• On Zen 3 and Zen 4 it can go from 13 to 7. Nobody has timed those.
• On Zen 5 the PR as posted took it from 26 to 14, and that is what my tables showed.

The chain columns count what a step must wait for from the step before. The fold still runs 12 maxes in a row inside a step, but the next step does not wait for them.

These counts are my arithmetic from the first table, and they are floors for the vector work alone; the loop has other costs.

So the tables in the description are Zen 5 numbers. I should have named the machine, and tried a second compiler, before posting them.

Sources: the uops.info pages for VPMAXUB, VPSUBSB, PMAXUB and PSUBSB; the InstLatx64 dumps for EPYC 7742 and Threadripper 3970X (Zen 2), Ryzen 5000 (Zen 3), Ryzen 7000 and EPYC 9654 (Zen 4), Ryzen 9 9950X (Zen 5) and the Alder Lake P-core.

Why long models got slower

The widest kernel holds 14 band vectors, xEv and beginv, which is all 16 vector registers. The second maximum was a 17th value. gcc 15 puts two values on the stack in the loop that short sequences run in. gcc 13 puts all 16 there: 31 loads and stores per step where main has none. I can reproduce that here with gcc 13.4. For the 3,841-position model at L = 30, the gain falls from 46% to 24%. On a CPU that gets no gain from the change, that is a net loss, and your slow cells are the ones where that loop dominates.

What I've changed

I've pushed two commits in place of the old one.

1. Fold. Each step folds its vectors into a temporary and joins that to xEv once, so only one max per step waits on the step before. It is the same number of operations as main, and there is no second long-lived register. With gcc 13.4 and gcc 15.2, the loops that step a band have the same stack traffic as main at every width but 14, where one value stays on the stack. A band of two vectors keeps the plain chain.
2. A gcc pragma. While counting registers I found that gcc, with -ftree-coalesce-vars on (the default), copies every band vector to a scratch register and back at each step. In the SSE build, the loop for a 14-vector band is 65 instructions per step, 28 of them register copies; without the copies it is under 40. main has this too, and the scratch register is what leaves the wide bands short of registers. #pragma GCC optimize ("no-tree-coalesce-vars") in the three kernel files removes the copies and most of the stack traffic, including that of main's long-sequence loop. It is a separate commit.

In short, for models of 100 positions and up, as change in time per call against main:

CPU	pragma	pragma with the fold	basis
Intel	no change, or a small gain	the same: the fold adds nothing there	a model for the pragma; unit counts above for the fold
Zen 2	up to 35% faster; 7 of 144 cells slower, worst by 9%	up to 40% faster; 1 of 144 cells slower, by 6%	measured, below
Zen 3, Zen 4	not known	should gain: same units as Zen 5	not measured; unit counts above for the fold
Zen 5	up to 21% faster; 7 of 144 cells slower, worst by 3%	20% to 72% faster in every cell	measured, below

The measured figures cover both compilers and all three kernels, 144 cells per CPU. A cell counts as slower from 1.5%; smaller differences are within noise.

On my machine, with gcc 13.4. Each cell is the pragma alone on main, then the fold, then the fold with the pragma, as change in time per call against main. gcc 15.2 gives the same picture.

Zen 5, AVX2:

model	30	100	300	1000
100	−2% / −18% / −20%	−1% / −33% / −34%	0% / −37% / −37%	0% / −39% / −39%
200	−5% / −39% / −42%	−2% / −55% / −56%	−1% / −59% / −59%	0% / −60% / −61%
400	−2% / −49% / −50%	−2% / −59% / −60%	−3% / −62% / −64%	−3% / −63% / −65%
1000	−2% / −48% / −57%	−16% / −55% / −68%	−19% / −57% / −71%	−21% / −57% / −72%
2000	0% / −54% / −59%	−7% / −47% / −63%	−17% / −49% / −69%	−20% / −50% / −71%
3841	+3% / −57% / −60%	+2% / −54% / −61%	−13% / −48% / −66%	−20% / −47% / −69%

Zen 5, SSE-only build:

model	30	100	300	1000
100	−9% / −42% / −47%	−3% / −58% / −60%	−1% / −62% / −63%	0% / −63% / −64%
200	−1% / −48% / −51%	+1% / −59% / −62%	+3% / −61% / −65%	+3% / −63% / −66%
400	−6% / −46% / −56%	−2% / −48% / −62%	0% / −48% / −64%	0% / −48% / −65%
1000	−2% / −53% / −60%	0% / −44% / −61%	0% / −40% / −63%	−1% / −39% / −64%
2000	−1% / −56% / −63%	+1% / −58% / −64%	0% / −44% / −64%	−1% / −37% / −65%
3841	0% / −56% / −61%	−1% / −54% / −61%	−1% / −46% / −60%	0% / −36% / −61%

I also ran all of this on a Zen 2 machine, a Threadripper 3990X, with gcc 13.4 and gcc 15.2: the same design as your benchmark, one pinned core, change against main. With gcc 13.4 the PR as first posted shows your slowdown there: +26% for the 1,000-position model at L = 30 (you had +26%) and +47% for 3,841 (you had +37%), in the AVX2 kernel. So I think the two benchmarks are comparable.

Zen 2 with gcc 13.4. Each cell is the pragma alone on main, then the fold, then the fold with the pragma.

Zen 2, AVX2:

model	30	100	300	1000
100	0% / −1% / −2%	0% / 0% / 0%	0% / 0% / 0%	0% / 0% / 0%
200	−4% / −3% / −7%	−2% / −4% / −9%	−1% / −4% / −7%	−1% / −4% / −9%
400	−3% / −6% / −11%	−4% / −10% / −15%	−4% / −15% / −16%	−5% / −15% / −16%
1000	−16% / +4% / −16%	−29% / −12% / −32%	−34% / −20% / −38%	−35% / −22% / −40%
2000	−3% / 0% / −9%	−20% / −1% / −23%	−32% / −6% / −35%	−35% / −7% / −38%
3841	−1% / 0% / −8%	−4% / +1% / −8%	−25% / +1% / −26%	−32% / 0% / −32%

Zen 2, SSE-only build:

model	30	100	300	1000
100	−5% / −3% / −7%	−2% / −4% / −11%	0% / −3% / −11%	−1% / −4% / −12%
200	−8% / −3% / +1%	0% / −9% / −14%	+3% / −10% / −19%	+6% / −11% / −20%
400	−9% / +6% / −4%	−12% / +3% / −16%	−12% / +2% / −22%	−12% / +1% / −24%
1000	−11% / +6% / −10%	−19% / +1% / −18%	−21% / +1% / −28%	−22% / +2% / −30%
2000	−12% / +2% / −13%	−8% / +2% / −15%	−22% / +2% / −26%	−27% / +1% / −32%
3841	−9% / +1% / −10%	−9% / +3% / −10%	−18% / +2% / −19%	−29% / +1% / −30%

So on Zen 2 the pragma does most of the work, and the fold by itself does little. The two together are best in the AVX2 kernel and on long sequences. On short sequences in the SSE-only build the pragma alone is as good or slightly better. Across both compilers and all three kernels in this run, the slowest cell with both changes is +6% (model 200 at L = 30, in the SSE kernel of the AVX2 build). The pragma alone has a cell at +9%, and the fold alone has one at +13%.

One caution: single cells on that machine are not stable from run to run. One cell of the fold with the pragma read +88% in one run and +1% in the next, with rounds inside each run agreeing to 1%; I could not find the cause, so a stray slow cell is worth a second run.

What I checked

• Return status and score match main on 22,000 calls per build, for the AVX2 build and the SSE-only build, with gcc 13.4 and gcc 15.2, on both machines.
• bathsearch hit tables match main for one model of 3,841 positions and for 50 models, in both builds.
• make check passes on both builds.

If you have time to run your benchmark on the new head, I'd like to know whether the slowdown is gone on your machines. I'll update the title and description once you've had a look.

🤖 Generated with Claude Code

[kevinmoran1988]

Problem. The SSV filter spends nearly all its time in band kernels. Each step of a band does one saturating subtract and one byte max per vector, 26 operations for a band of 13, and in main every max needs the result of the max before it, with the first max of a step needing the last max of the step before. So a step cannot start until the one before has finished all 13. PR 36 first split that chain in two. On my Zen 5 it was 35–47% faster, but you measured no gain on Intel, little on Zen 2, and a 13–61% slowdown for long models on short sequences. The slowdown had its own cause: the second running maximum was a 17th value for 16 vector registers, and gcc 13 responded by pushing every band vector to the stack at each step.

Why it depends on the CPU. Two things limit a step: the chain, and how many execution units the core has for these operations. Intel runs both operations on the same two units, so 26 operations take 13 cycles, exactly as long as the chain. It was never waiting on the chain, and shortening it gains nothing. AMD's cores have more units than the chain let them use. Zen 3, 4 and 5 run byte max on four units and saturating subtract on two, enough for about 7 cycles a step, but the chain held them to 13. Zen 5 was hit hardest, because a chained vector operation costs it 2 cycles where every other core pays 1: a 26-cycle chain for work its units could do in 7.

Solution. Two changes replace the original. The fold combines a step's vectors into a short-lived temporary and joins the running maximum once. The maxes inside a step still follow one another, but only that one join waits on the step before, so steps overlap, and no extra long-lived register is needed.

main:   step 1 [=============]
        step 2                [=============]
        step 3                               [=============]

fold:   step 1 [============]x
        step 2        [============]x
        step 3               [============]x

x is the one max that joins the running maximum. On Zen 3 to Zen 5 a new step can start every 7 cycles. On Intel the units only allow one every 13, so nothing changes there. A step then takes what the units allow: 26 to 7 on Zen 5, 13 to 7 on Zen 3 and 4, 13 to 9 on Zen 2, and 13 to 13 on Intel.

Separately, a gcc pragma turns off variable coalescing in the three kernel files, because gcc otherwise copies every band vector through a scratch register at each step (28 wasted copies out of 65 instructions in the widest band) and loses a register doing it; that saving does not depend on units. Together they measure 20–72% faster on Zen 5 and up to 40% faster on Zen 2, where the pragma does most of the work and the worst cell is +6%. Intel should see no change from the fold and is unmeasured for the pragma, and Zen 3 and Zen 4 are predicted from their unit counts, not timed. Scores match main in every test, and make check passes.

[kevinmoran1988]

I've added a third commit: the fold now applies with gcc only.

I timed clang, which I had not done before (clang 21, on the Zen 5 machine). Two things came out of it.

clang did not need the fold. It reorders a step's maxes by itself and has no scratch copies (what the gcc pragma fixes), so main built with clang is already about as fast as gcc with both commits. Time per call in the AVX2 kernel:

model × sequence	gcc 13.4, main	gcc 13.4, both commits	clang 21, main
200 × 300	734 ns	300 ns	316 ns
1000 × 1000	14,079 ns	4,024 ns	4,234 ns
3841 × 1000	53,692 ns	16,826 ns	16,779 ns

The fold made clang slower in places. Against clang-built main, over the 24 cells for models of 100 positions and up, 15 were slower in the SSE kernels, by up to 22%, and 7 were slower in the AVX2 kernel, by up to 8%. The slower cells are mostly those whose bands are 13 or 14 vectors wide. Narrower bands were faster, by up to 18%. I'm looking into why some bands speed up while others do not, but for now the fold is guarded. The pragma was never a problem there: clang skips it.

So both changes are now fixes for how gcc compiles this kernel. With clang 21 the three files compile to the same machine code as main. With gcc 13.4 and gcc 15.2 they compile to the same machine code as before this commit, so the gcc numbers above stand. Scores still match main, and make check passes with gcc and with clang, for the AVX2 build and the SSE-only build.

I have only timed clang on Zen 5, and only that one version.

Update: older gcc. Many clusters run older versions of gcc, so I tested more of them: 8.5, 11.5, 12.5 and 14.3, on top of 13.4 and 15.2. gcc 9 and 10 are still untested.

• gcc 11 to 15, on the Zen 5 machine: the commits are faster than main in every cell of the grid (models of 100 positions and up, all three kernels) and for every band width from 2 to 14. The five versions agree to within about 3 points.
• gcc 8.5, on a cluster node with an EPYC 9654 (Zen 4): the fold is 9–41% faster in the AVX2 kernel and up to 45% faster in the SSE kernels, with no cell more than 1% slower. The pragma is another matter there. By itself it was up to 23% slower than main, and on top of the fold up to 8% slower.

So I've added a fourth commit: the pragma now applies from gcc 11 on, and older gcc gets the fold alone. For gcc 11 to 15 the machine code is unchanged.

Scores match main in every one of these builds, on 22,000 calls each.

🤖 Generated with Claude Code