Skip to content

[auto-optimize]: Speed up mem seeding, SA lookup, extension and CIGARs - #656

Open
jmschrei wants to merge 1 commit into
fg-labs:mainfrom
jmschrei:auto-optimize/2026-10-27-port-A69
Open

jmschrei wants to merge 1 commit into
fg-labs:mainfrom
jmschrei:auto-optimize/2026-10-27-port-A69

Conversation

@jmschrei

Copy link
Copy Markdown

This optimization was found with auto-optimize, a tool developed by the Programmable Genomics Laboratory at UMass Chan Medical School: it makes bwa-mem3 mem faster in seeding, SA lookup, seed extension, mate rescue, reference fetch, CIGAR generation and start-up, with byte-identical output (see Timing for what was measured and on which base, and How the output was checked).

Changes

The changes are grouped by area, and within each group they are listed in the order they were made. Paths are relative to the repository root. No CLI option, CLI default, index format or build flag changes. One environment variable is new (BWA3_BWTSEED_JUMP), and three existing rescue switches now also select paths in the AVX-512BW rescue kernel.

Seeding: the SMEM and bwtSeed passes

  • Branch-free BWT symbol (src/FMI_search.cpp, FMI_search::call_one_step). The BWT symbol at sp is computed without branches.
  • Fewer look-ahead prefetches (src/FMI_search.cpp, src/FMI_search.h, the lockstep SMEM and bwtSeed drivers). The SMEM driver drops its cross-slot T1 look-ahead and the bwtSeed driver its T1 look-ahead.
  • Backward step in registers (src/FMI_search.cpp, src/FMI_search.h, ls_advance_backward_step, FMI_search::backwardExt_konly_cp). The backward step keeps its slot's fields in locals, and its loops extend through backwardExt_konly_cp, which takes cp_occ and count[a] as arguments, so they stay in registers across the prev[] stores.
  • Forward phase writes prev[] downward (src/FMI_search.cpp). The forward phase fills prev[] from the top instead of filling it upward and reversing it.
  • Unique-match walks compare against the text (src/FMI_search.cpp, src/FMI_search.h). A walk whose interval is down to one occurrence (s == 1) with 16 or more bases left compares the read with the reference text 8 bases at a time instead of stepping the FM index, and a walk's last bases at s == 1 take a fast path.
  • The bwtSeed pass reuses the SMEM pass's first walk (src/FMI_search.cpp, src/FMI_search.h, src/bwamem.cpp, BsdFirstWalk, bwtSeedStrategyAllPosOneThread_lockstep, mem_collect_smem). A read's first bwtSeed walk (from base 0) is a prefix of the SMEM pass's first forward walk. The SMEM walk records where the bwtSeed walk would stop and its seed, and the bwtSeed pass emits that seed and continues from there instead of walking it again.
  • The SMEM driver steps its slots by kind (src/FMI_search.cpp, getSMEMsOnePosOneThread_lockstep). Slots sit on forward, unique-match and backward lists and are stepped list by list instead of switching on each slot's phase. The driver reads LockstepParked::cursor().
  • Interval table for the first bases of a walk (src/FMI_search.cpp, src/FMI_search.h, load_index, load_index_from_shm). At index load, the bidirectional interval of every string of 1 to 10 bases is tabulated (33.6 MB, about 2 ms with the load's threads). Each lockstep bwtSeed walk starts at the table's interval for its first bases, and each lockstep SMEM forward walk reads its first intervals from the table in one visit, testing each as a step would. BWA3_BWTSEED_JUMP=<d> sets the depth (2 to 12; 0 for no table).
  • bwtSeed walks that cannot reach -k are retired (src/FMI_search.cpp, src/FMI_search.h). A bwtSeed walk that can no longer reach the minimum seed length stops early.

SA lookup

  • 64 SA-resolve lanes on x86 (src/FMI_search.cpp, SA_RESOLVE_LANES_DEFAULT). The x86 default goes from 40 to 64. The comment records the sweep it came from.
  • Unconditional prefetches (src/FMI_search.cpp, FMI_search::get_sa_entries_prefetch). Both prefetches are issued without a test.
  • Single-occurrence SMEMs keep their text position (src/FMI_search.cpp, src/FMI_search.h, src/bwamem.cpp, mem_collect_smem). An SMEM with one occurrence keeps the text position its LF walk resolved, so SA resolution writes it without a second walk. A contained single-occurrence SMEM takes the position of the SMEM that contains it. The tagging runs over each read's SMEMs after smem_sort_rid_mn and before smem_dedup.

Seed extension

  • AVX-512BW 8-bit lean row (src/bandedSWA.cpp, src/bandedSWA_compact512.inc, smithWaterman512_8Impl and its lane-compaction driver). The AVX-512BW 8-bit kernel and its compaction driver use the aarch64 kernel's row: a counted band trim, a narrow per-row epilogue with a gated z-drop, band narrowing from offsets that the cell loops record, and 64-lane vector reductions. It works in both score domains of the kernel (perf(bsw): keep the AVX-512BW 8-bit extension in a biased signed domain #647): its zero tests, the z-drop gate and the result bytes read the score through bsw8_hf_or512 / bsw8_unbias512. The z-drop gate counts live lanes.
  • Probe band for longer non-HIT pairs (src/bandedSWA.h, src/ungapped_ext.h, src/bwamem.cpp, SeqPair::probe_band, ungapped_analyze_pb, bsw_tb_probe_rung). ungapped_analyze_pb records a probe band in the new SeqPair::probe_band field, so the certified probe rung also takes non-HIT pairs with a query of 129 to 512 bases.
  • Probe at clip slack (the same files and functions). A non-HIT pair whose band numerator is not positive is probed at its clip slack.
  • 8-bit re-sort keyed on the query length (src/bwamem.cpp, sortPairsLen, PairSortKey). Before the 8-bit kernel, pairs are sorted by len2 instead of max(len1, len2), since the 8-bit kernels bound each lane by its query; a call diverted to the 16-bit kernel keeps max(len1, len2).
  • Stable tier partition (src/bwamem.cpp, sortPairsLenExt). The extension pairs are split into the 8-bit, 16-bit and scalar tiers by a stable partition, with the 8-bit pairs compacted in place, instead of a counting sort by minval (three passes over the pairs, a scatter and a copy back, and a 1,088-bin histogram on every call). The kernels' own sorts decide the packing order, so only their tie-break changes.
  • Ungapped tail in one step (src/ungapped_ext.h, ungapped_analyze_pb). The last N mod 16 bases are handled in one 16-byte AVX2 step.
  • Mismatch walk over set bits (src/ungapped_ext.h, src/bwamem.cpp, ungapped_walk_mis). The walk iterates the set bits of each mismatch word.
  • Fewer, fuller probe rung buckets (src/bwamem.cpp, bsw_tb_probe_rung). The certified probe rung's bucket ceilings go from {4, 8, 12, 16, 20, 24, 28, 30} to {8, 16, 24, 30}, and below the top admissible ceiling each bucket's partial last lane group joins the next ceiling, so fewer kernel calls run, each with fuller groups.

Mate rescue

  • Two-row sweep in the AVX-512BW 8-bit rescue kernel (src/kswv.cpp, src/kswv.h, src/rescue_env.h, kswv512_u8_impl). The kernel's FScan body sweeps two target rows per pass, as the NEON kernel does, sharing column loads and the vertical carry. BWA3_RESCUE_ROWPAIR=0 selects one row at a time.
  • All-pad columns and a reverse query-end rescan (src/kswv.cpp, src/rescue_env.h, kswv512_u8_impl). In the signed-domain FScan bodies, the query columns where every lane holds the pad code (X86_QPAD8) skip the score gather, since the pad's score is 0; the query-end rescan scans a block last to first. BWA3_RESCUE_PADCOLS=0 and BWA3_RESCUE_QE_REVSCAN=0 select the previous forms.

Chaining, reference windows and CIGARs

  • Chain windows hinted 4 chains ahead (src/bwamem.cpp, mem_chain2aln_across_reads_V2). The extension driver hints each chain's reference window 4 chains ahead, in flat read-then-chain order across reads, instead of one chain ahead.
  • mem_chain_flt sorts (w, index) pairs (src/bwamem.cpp, mem_chain_flt, chain_flt_sort). From 32 chains up, the sort by weight runs on 8-byte (w, index) pairs with the same algorithm and comparator, which gives the same permutation, ties included, and the 48-byte chains then move once.
  • 32-byte reverse copies (src/ext_reverse_copy.h, ext_reverse_copy, ext_reverse_inplace_avx2). The reverse copy runs in 32-byte AVX2 blocks.
  • No per-seed log in the flat chaining index (src/bwamem.cpp, ChainFlat, chain_add_seed, chain_replay_read_prefix). The flat index (perf(chain): pack the flat chaining index into one (pos, chain index) key array #642's packed keys) keeps no copy of a read's seeds. When it declines a read, the default path resolves the read's seeds again up to the declining one and chains them through the B-tree; the opt-in seed reorder modes replay their already ordered seeds.
  • Reference windows written in 16-base pieces (src/bntseq.cpp, bns_get_seq_into). The start and end of a window are written as 16-base pieces at any base phase, on both strands.
  • CIGAR mismatches 32 bytes at a time (src/bwa.cpp, bwa_gen_cigar4). An M op's mismatches are found 32 bytes at a time.
  • In-place AVX2 reversals in bwa_gen_cigar4 (src/bwa.cpp, src/ext_reverse_copy.h, bwa_gen_cigar4). The query and reference windows are reversed in place with AVX2.
  • Chain window from the look-ahead (src/bwamem.cpp, mem_chain2aln_across_reads_V2). The look-ahead computes each chain's union window over all its seeds and its largest n_hits, hints that window, and passes both to the window build through a ring, so the build does not walk the seeds again.

Start-up

  • Inputs read while the index loads (src/fastmap.cpp, src/fast_reader_bseq.c, src/fast_reader_bseq.h, main_mem, mem_open_inputs, preread_main, bseq_diag_t). main_mem opens the inputs and reads the first slice on a helper thread while the index loads. The reader's warnings and a malformed-input error are held back and replayed after the load, so what a run prints and its exit status stay as before, apart from what Behaviour changes lists. Open failures, including fix: report corrupt index, bad -T and stream/pthread failures as errors #649's --legacy-reader gzdopen checks, are reported after the load as before, and a failed join of the helper thread exits as fix: report corrupt index, bad -T and stream/pthread failures as errors #649's pipeline joins do. Under --legacy-reader the thread only opens the inputs.
  • Tick rate measured over the run (src/main.cpp, src/profiling.cpp, src/profiling.h, proc_freq_mark_start, proc_freq_resolve). The __rdtsc() tick rate is measured between start-up and the first report that needs it, instead of over a 5 ms sleep at start-up.

Documentation

  • Knob tables (src/rescue_env.h, docs/_generated/rescue/knobs.md, docs/src/developer-guide/rescue-kswv.md, docs/src/whats-different/performance.md). They list BWA3_BWTSEED_JUMP and the AVX-512BW paths of BWA3_RESCUE_ROWPAIR, BWA3_RESCUE_PADCOLS and BWA3_RESCUE_QE_REVSCAN.
  • Flat chaining replay (docs/src/whats-different/equivalence.md). The flat chaining section describes the replay without the seed log.

Tests

The new unit-test suite is in #654. All runs below used clang 19 and OMP_NUM_THREADS=4, each in a fresh export of the tree.

run tree unit cases other doctest binaries mem CLI scripts (checks) exit
make test main b5a4449 488/488 1/1, 19/19, 4/4 not in main 0
make test this branch 488/488 1/1, 19/19, 4/4 not in this branch 0
make arch=avx2 USE_MIMALLOC=1 test this branch 522/522 1/1, 19/19, 4/4 not in this branch 0
make test this branch merged with the suite 588/588 1/1, 19/19, 4/4 112, 315, 161, 165, 125 0
make arch=avx2 USE_MIMALLOC=1 test this branch merged with the suite 622/622 1/1, 19/19, 4/4 112, 315, 161, 165, 125 0

No test failed. The only skips are shell_lint and its selftest, because shellcheck and shfmt are not installed on the test host. On this branch and on main the unit binary makes the same 5,477,949 assertions, and upstream's test_lockstep_width.cpp subcase "shipped configuration: Zen 5 is 32/32", run alone, makes its 2 assertions and passes.

CI's AVX-512BW tier step, on a unit binary built at arch=avx512bw, passed: -ts=unit/kswv 38/38, the lane-compaction tests with BWA_REQUIRE_COMPACTION_TIER=avx512bw 9/9, the 8-bit z-drop gate tests 1/1, -ts=unit/bandedswa-signed8 with BWA_REQUIRE_SIGNED8=1 2/2, and the whole unit binary 522/522. These CI scripts outside make test also pass on this branch and on main: rescue_docs_lint.sh (26 knobs, generated table current) and its selftest, arg_range_validation.sh, chain_flat_replay_identity.sh in fixture mode (200 tie and 200 cap declines), rescue_prune_identity.sh, dash_read_as_n.sh and seedcov_mapq_legacy.sh. The chr22 jobs, the debug-macro, ASAN, MSan and coverage builds, and the ARM64 and macOS rows were not run.

gcc 15 builds this branch and main with the same two warnings ('_get_pac' redefined in src/bwamem.cpp, and -Walloc-size-larger-than= from src/kvec.h inlined into mem_chain_flt). With clang 19, both give the same 9 distinct warnings.

Behaviour changes

No SAM, BAM or index file changes (How the output was checked, below). What a user can see differ, in that comparison, is on stderr or in how some failures happen:

  • A new environment variable, BWA3_BWTSEED_JUMP, sets the depth of the interval table (2 to 12; 0 for no table). An invalid value prints an error line and is ignored.
  • The inputs are opened while the index loads. Messages written while opening them now appear during the index load instead of after it: a <cmd| input's shell errors, or connect: Connection refused for an unreachable http:// or ftp:// input. Under --legacy-reader, zlib's message for a corrupt gzip file names <fd:5> instead of <fd:4>. The exit status and the final error message are unchanged.
  • The interval table is allocated at index load (33.6 MB at the default depth). Under ulimit -v 300MB, a run now fails at index load with Allocation of 0.03 GB for seed_jump failed instead of later in worker_alloc. Both exit with an error.
  • An assertion's reported line moved. On a damaged index (a yeast BWT with E. coli .pac/.ann/.amb), both builds abort on the same assert in bns_clamp_window, reported at the assert's own line, which is bntseq.cpp:769 in main b5a4449 and :865 here.
  • BWAMEM3_DEDUP_STATS=1 counters differ, for example total_jobs=1425 distinct_jobs=1117 against main's 1423 / 1105 on one input, because the change generates a slightly different set of extension jobs internally. The SAM is identical.
  • "Shringking not supported yet." was printed once more here than in main on two inputs, steadily over five runs of each. Main itself prints this line a varying number of times from run to run ([auto-optimize]: Bugs found while testing bwa-mem3 b5a4449 #655, bug 6).

Timing

These timings were measured on the port of this change onto 6c2cd79, with the Zen 5 bwtSeed lockstep width at 48. This pull request does not include that width change, and the timings were not re-measured on b5a4449. The search's own sweep of the width found that 48 instead of 32 saved 0.48 of about 57 seeding thread-seconds over the three timed inputs, about 0.02 s on the mean below.

Setup:

  • Host: an AMD EPYC 9575F (Zen 5, AVX-512BW tier), so the width of 48 was in effect.
  • Builds: the Makefile's default build with clang 19.1.7, in a fresh clone; the base is main 6c2cd79.
  • Index: GRCh38 no-alt, in the page cache.
  • Runs: each is one bwa-mem3 mem -t 8 -K 80000000 process pinned to one NUMA node, timed from spawn to exit (index load, mapping, SAM written to a file). There were three runs per build, each on a different slice of every input, and the medians are reported.

Inputs:

  • WGS: 1.2M GIAB HG002 Illumina 2×148 pairs.
  • WES: 3M HG002 exome 2×151 pairs (SRA SRR14724472).
  • ChIP-seq: 9M ENCODE SE50 reads (ENCFF459LYK).
input threads main 6c2cd79 (s) this change on 6c2cd79 (s) difference
WGS 8 9.285 7.352 −1.932 s (−20.8%)
WES 8 8.794 6.170 −2.624 s (−29.8%)
ChIP-seq 8 8.826 6.685 −2.141 s (−24.3%)
mean of the three 8 8.968 6.728 −2.240 s (−25.0%, 1.333x)
WGS 16 4.983 3.956 −1.027 s (−20.6%)

Median CPU seconds fell on all three inputs at 8 threads: WGS 74.9 to 59.9, WES 72.4 to 51.8, and ChIP-seq 72.0 to 55.4. Peak RSS stayed within 0.11 GB of main on every case; the largest difference was ChIP-seq at 14.54 against 14.42 GB.

How the output was checked

The timings above measure speed. Separately, we ran main and this change side by side on a collection of public sequencing data and hand-made malformed files. That comparison used the port of this change onto 6c2cd79, with the bwtSeed width at 48, against main 6c2cd79; this commit itself was checked on the search's benchmark cases (This commit, below).

  • Scale: 16,991 cases (a command line, its environment and its input), 16,059 of them run again with both builds forced to the AVX2 tier. That is 33,050 comparisons: 72,180 process runs and 27.1 hours of process time, including the controls.
  • What was compared: for every case, the exit status, stdout, stderr and every file written, byte for byte. Both builds were run by the same name from the same directory, so @PG lines are compared too. Before stderr was compared, timing lines, thread progress lines and the lockstep width line (32 in main, 48 in the build compared) were replaced by tokens.
  • Threads: at a fixed -K, bwa-mem3's output does not depend on the thread count, so every case ran at -K 10000000 (the full-depth files at -K 80000000), except the cases that test the default, -K and --chunk-cap.
  • Inputs: 1,439 files: reads from 53 assay directories (human, mouse, yeast, fly, Arabidopsis, bacterial mock communities and SARS-CoV-2), the collection's malformed FASTQ and FASTA sets, and 35 genomes and reference sequences. Human reads, including those released against hg19 and hs1, were aligned to GRCh38 (no-alt). Other species were aligned to small indexes that both builds built, and those indexes were compared too.

Builds compared

Pair main 6c2cd79 This change on 6c2cd79 Why
Native make make what users run; on this CPU both run their AVX-512BW tier
Forced AVX2 the same binary with BWAMEM3_FORCE_TIER=avx2 the same the AVX2 kernels, which this CPU otherwise never selects; the change rewrote SIMD code
Controls make main again whether a difference also appears between two runs of main: 3,313 cases sampled across the suites, and 905 cases rerun after they differed

Both builds are clang 19.1.7 builds of the Makefile's default goal, the multi-tier binary at -O3 -flto=full.

What ran

Suite Cases (native / AVX2) What
Read tiers 753 / 753 every read file of every input class in the collection at tiny and small size, and one per class at medium and large: short reads, 10x, proximity ligation (Hi-C, Micro-C, Capture-C, Pore-C), bisulfite, long reads, interleaved pairs and FASTA queries; mem -K 10000000 at 1 to 4 threads, -p on interleaved files
Options 14,806 / 14,692 every row of the command-line surface of mem, index, re-sa, version and the top level (556 rows: each option and each environment variable mem reads, with several values) on each accepted input class's small file
Combinations 137 / 70 a pairwise covering array over 13 factors (output format, mode including --fast and --compat, pairing, split reads, XA, scoring, seeding, dedup, reader, batch size, SIMD tier, header and threads) on medium files: every pair of values of two factors occurs in some case
Threads 32 / 32 1, 2, 3, 4, 7, 8 and 16 threads, the same count with and without -K, and this change at N threads against main at 1
Index builds 268 / 22 index on 35 genomes and reference sequences from SARS-CoV-2 and phage lambda to fission yeast, chr22, chr19 and chrM; each index option; re-sa; the two --meth indexes; mem on each build's own index
Full depth 22 / 22 22 full-size read files (9.0 GB of FASTQ and FASTA) from 15 assays, including Micro-C, NOMe-seq, MPRA, perturb-seq, ancient DNA, amplicon long reads and hs1 ChIP-seq, at -t 8 -K 80000000
Edge-case files 685 / 242 602 edge-case files: 155 FASTQ and FASTA query files through mem (truncated or empty files, CRLF line ends, mates out of sync or of unequal counts, gzip CRC errors, N-only barcodes, and others), 447 malformed FASTA files through index
Invalid input 207 / 167 out-of-range and malformed values of every validated option, refused combinations, unknown options, missing, empty, unreadable and directory paths, BAM and CRAM given as reads, bwa 0.7 indexes, http:// and ftp:// inputs that cannot connect, bad environment values
Damaged indexes 58 / 40 copies of a sacCer3 index with .bwt.2bit.64 truncated at 0, 8, 48 and 4,096 bytes, at half and at one byte short, a flipped bit, missing, empty or garbage .pac/.ann/.amb, files from another genome, and a damaged --meth seed index; mem and re-sa on each
Faults 21 / 19 SIGTERM, SIGINT and SIGHUP mid-run, SIGPIPE, read-only and missing output directories, ulimit -f and ulimit -v, closed stdout, reads on stdin, FIFOs, `<cmd
shm 2 / 0 shm help only (What this testing does not cover)

Results

Pair Cases Identical Same failure Differ
Native 16,991 15,380 1,111 500
Forced AVX2 16,059 14,968 606 485
Controls (main against main) 4,218 3,115 188 915

"Identical" is valid input with the same exit status, stdout, files and stderr. "Same failure" is malformed input, or two runs that both failed, with the same exit status, files and stderr. No SAM, BAM or index file differs between the builds in any case, valid or malformed, on either tier. Every difference is on stderr, in how a malformed or fault case fails, or in a debug file that several threads write. Each cause, over both pairs, with the evidence for it:

Cause Cases Evidence
"Shringking not supported yet." printed a different number of times 888 main prints it a varying number of times between runs (#655, bug 6); on 2 cases the count was steady over five runs of each build and one higher here (Behaviour changes)
BWAMEM3_DEDUP_STATS=1 counters 26 22 vary between runs of main itself; 4 were steady in each build and different (Behaviour changes)
-v 4 and -v 5 debug lines from worker threads, interleaved mid-line 20 main differs from itself the same way (#655, bug 7)
reader progress lines and timing probes in another order or with other values 10 main differs from itself the same way in 6; the other 4 are two progress lines in another order on a file with unequal mate counts, with the same exit status and final error
line order of the BWAMEM3_DUMP_PAIRS debug file 10 written by several threads; main differs from itself in 9, and the 10th matched main in 4 of 5 runs
BWA3_BWTSEED_JUMP set to an invalid value 10 the variable is new; main does not read it
which misnamed or unpaired read a multithreaded run reports first 6 main differs from itself the same way in 5; the 6th, on a file with mates out of sync, has the same exit status and final error
messages from opening the inputs printed earlier (a `<cmd ` input's shell errors, an unreachable URL) 6
an assertion's source line number 4 the same assert, moved by lines added above it
zlib's file descriptor number under --legacy-reader 2 the inputs are opened during the index load
ulimit -v 300MB fails at index load 2 the interval table is allocated there; both builds exit with an error
-c 2147483647 1 both builds segfault on it, main in 4 of 5 runs and this change in 5 of 5 (#655, bug 2)

On the two cases where both builds failed on valid input (counted under "Same failure"), a Pore-C file at -K 10000000 -t 4 stops both at the read-buffer capacity check, which tells the user to lower -K.

Memory

Peak resident memory on the GRCh38 index is unchanged in practice. Over the 11,393 identical native cases that peaked above 5 GB, the ratio of this change to main ranges from 0.983 to 1.051, with a median of 1.003 (0.03 GB). The largest increases:

  • BWA3_BWTSEED_JUMP=12, the deepest interval table, which grows fourfold per base of depth: about 0.5 GB more, in 4 cases.
  • A full-depth MPRA file (PE15, -t 8): 16.2 GB against 15.6 GB in main.
  • PacBio HG002 long reads at the small tier: 0.41 GB more on 13.1 GB.

On smaller indexes, the 246 identical cases that peaked between 0.5 and 5 GB range from 0.958 to 1.225, at most 0.12 GB more.

This commit

The commit in this pull request was run with main b5a4449 on the search's own benchmark cases and compared under the search's rule: SAM byte-identical except the @PG lines, BAM compared decompressed without @PG, and every other file byte for byte. There were 59 case-and-input-slice pairs: 8 cases at 8 or 16 threads on three slices each, and 35 option cases on one. These cover --fast, both --compat targets, -M, -Y, -a, -C -R -H, -p, --bam and --bam=6, --hic, -I -S -P, scoring and seeding options, the dedup modes, --legacy-reader, SAM on stdout, small -K, five cases forced to the AVX2 tier, and index, re-sa and --meth on a small genome. All 59 match, and every process exited 0. Main's outputs on these cases did not change from 6c2cd79 to b5a4449.

What this testing does not cover

  • This commit in the large comparison: the comparison above ran on the 6c2cd79 port. This commit differs from it in the lockstep width, which does not change output, and in two pieces re-implemented on b5a4449's code: the flat chaining index without its seed log (perf(chain): pack the flat chaining index into one (pos, chain index) key array #642) and opening the inputs during the index load (fix: report corrupt index, bad -T and stream/pthread failures as errors #649's gzdopen checks). Those were checked by the 59 cases above and the tests, not by the large comparison.
  • Platforms: only x86-64 Linux, one AMD Zen 5 CPU and clang 19. We compared no Intel CPU, ARM (NEON), macOS or other compiler; gcc 15 was only built and tested (Tests).
  • Sanitizers: no sanitizer build was part of this comparison.
  • Not run: the shm suite, which stages indexes in /dev/shm; whole-genome indexes other than GRCh38 (mouse, zebrafish) and a GRCh38 --meth index; human reads against hg19 and hs1 indexes; Phred+64 paired-end reads; genomes above 4.3 Gbp, and the five barley edge files. Each option ran on one input size (small); the larger and other tiers of the collection were not crossed with every option.
  • Timing: this comparison measured no speed.

Bugs in main found along the way

We reproduced 7 bugs in main and filed them as #655; none comes from this change. Among them:

  • mem -u writes its widened alternative hits under XA:Z, where the help, the CLI docs and bwa write XB:Z.
  • mem -c 2147483647 segfaults.
  • mem -k 0 aborts on a pair whose mate is all N.
  • The 16-bit banded Smith-Waterman kernels read an unaligned stack array with an aligned SIMD load, which crashes in a gcc 15 -O0 build.

Port

This is a port of the search's measured winner onto b5a4449. The search ran on 0.15.1 (51ee7eb). The change was then ported onto 6c2cd79, where it was timed and validated as above, and now onto b5a4449. Where upstream had rewritten code the change touched, the change was re-implemented on upstream's code:

Mate rescue's in-place reverse keeps #644's seq_reverse_inplace. The search's AVX2 version of it was dropped.

The search also raised the Zen 5 bwtSeed lockstep width (BWA3_ZEN5_BWTSEED_LOCKSTEP_N) from 32 to 48. That change was left out, so the width stays 32 and src/lockstep_width.h is unchanged.

Related

🤖 Generated with Claude Code

The measured winner of a speed search on 0.15.1 (51ee7eb), ported onto
main b5a4449. The change keeps `bwa-mem3 mem` output byte for byte
(SAM, and BAM decompressed, without the @pg lines) and changes no CLI
option, CLI default, index format or build flag. The search also raised
the Zen 5 bwtSeed lockstep width from 32 to 48; that change is left out,
and BWA3_ZEN5_BWTSEED_LOCKSTEP_N stays 32.

SMEM seeding and the FM index (FMI_search.cpp/.h):
- call_one_step computes the BWT symbol at sp without branches.
- The SMEM lockstep driver steps its slots by kind (forward,
  unique-match and backward lists) instead of switching on each slot's
  phase; ls_advance_backward_step keeps the slot fields in locals; the
  forward phase writes prev[] downward instead of reversing it; a walk
  that reaches a unique match (s == 1) with 16 or more bases left
  compares the read with the text 8 bases at a time, with an s == 1
  fast path for a walk's last bases; and forward walks take their first
  steps from an interval table that load_index builds once per process
  (BWA3_BWTSEED_JUMP sets its depth, default 10, 0 for none).
- The bwtSeed pass starts each read where the SMEM pass's first walk
  recorded it would end, starts each walk at the table's interval of
  its first bases, and retires walks that cannot reach the minimum seed
  length.
- An s == 1 SMEM keeps the text position its LF walk resolved, so SA
  resolution writes it without a second walk; contained s == 1 SMEMs
  take it after the per-read sort (smem_sort_rid_mn).
- The SMEM and bwtSeed drivers drop two T1 look-ahead prefetches,
  get_sa_entries_prefetch issues its two prefetches unconditionally,
  and SA_RESOLVE_LANES_DEFAULT on x86 is 64.

Extension and mate rescue (bandedSWA.cpp/.h, bandedSWA_compact512.inc,
kswv.cpp/.h, rescue_env.h, ungapped_ext.h, bwamem.cpp):
- The AVX-512BW 8-bit extension kernel and its lane-compaction driver
  use the aarch64 kernel's lean row (counted band trim, narrow per-row
  epilogue with a gated z-drop, band narrowing from offsets recorded by
  the cell loops, vector reductions), in both score domains of the
  kernel: its zero tests, z-drop gate and result bytes read the score
  through bsw8_hf_or512 / bsw8_unbias512.
- kswv512_u8 runs the NEON kernel's two-row sweep, skips the score
  gather in all-pad query columns, and rescans query-end blocks last to
  first; BWA3_RESCUE_ROWPAIR, BWA3_RESCUE_PADCOLS and
  BWA3_RESCUE_QE_REVSCAN now select these paths there too.
- The 8-bit kernels' per-call re-sort keys on the query length,
  sortPairsLenExt partitions tiers stably, the certified probe rung uses
  ceilings {8, 16, 24, 30} and also takes non-HIT pairs of length
  129-512 and pairs whose band numerator is not positive.
- ungapped_analyze_pb handles the last N mod 16 bases in one 16-byte
  step and ungapped_walk_mis iterates the set bits of each mismatch
  word.
- mem_chain_flt sorts (w, index) pairs; ChainFlat no longer keeps a
  per-seed log, and a read it declines is replayed by resolving its
  seeds again.

Chains, reference fetch and CIGARs (bwamem.cpp, bntseq.cpp, bwa.cpp,
ext_reverse_copy.h):
- mem_chain2aln_across_reads_V2 hints each chain's reference window 4
  chains ahead and computes each chain's union window once.
- bns_get_seq_into writes the start and end of a window as 16-base
  pieces at any base phase, on both strands.
- bwa_gen_cigar4 finds an M op's mismatches 32 bytes at a time and
  reverses the query and reference window with AVX2; ext_reverse_copy
  copies in 32-byte AVX2 blocks.

Start-up (fastmap.cpp, fast_reader_bseq.c/.h, main.cpp,
profiling.cpp/.h):
- main_mem opens the inputs and reads the first slice on a helper
  thread while the index loads, replaying the reader's diagnostics in
  order.
- The __rdtsc() tick rate is measured over the run
  (proc_freq_mark_start / proc_freq_resolve, now in profiling.cpp)
  instead of over a 5 ms sleep at start-up.

Docs: the knob tables (rescue_env.h, docs/_generated/rescue/knobs.md,
performance.md, rescue-kswv.md) list BWA3_BWTSEED_JUMP and the
AVX-512BW paths of the three rescue knobs; equivalence.md describes the
flat chaining replay without the seed log.

Changes against the measured change, from what main changed since
51ee7eb:
- The 8-bit lean row is re-implemented on fg-labs#647's score-domain template
  (the zero tests, z-drop gate, gscore fold and result bytes above).
- The contained-SMEM tagging runs after fg-labs#648's smem_sort_rid_mn instead
  of inside the per-read sort loop fg-labs#648 replaced.
- Mate rescue's revseq keeps fg-labs#644's seq_reverse_inplace; the AVX2
  revseq of the measured change is dropped.
- The SMEM driver reads LockstepParked::cursor() (fg-labs#629), and the
  two-row kswv sweep tests fg-labs#636's X86_QPAD8 query pad.
- ChainFlat drops its seed log on fg-labs#642's packed (pos, chain index)
  keys.
- mem_open_inputs, where the inputs are now opened, carries fg-labs#649's
  checks of the --legacy-reader gzdopen calls, with the same messages
  and clean-up, and a failed join of the helper thread exits the way
  fg-labs#649 made the pipeline's joins exit.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@coderabbitai

coderabbitai Bot commented Oct 28, 2026 •

Copy link
Copy Markdown

Note

Reviews paused

Use the following commands to manage reviews:

  • @coderabbitai resume to resume automatic reviews.
  • @coderabbitai review to trigger a single review.

Use the checkboxes below for quick actions:

  • ▶️ Resume reviews
  • 🔍 Trigger review

Warning

Review limit reached

You've used all free OSS reviews for now. Wait for the free limit to reset to keep reviewing this public repository.

Next included review available in 48 minutes.

Check out review usage here.

View limit details

Limit details: You’ve used the included review currently available.

Learn how review limits work.

Review configuration:

Run configuration
Commits

Reviewing files that changed from the base of the PR and between b5a4449 and edb267b.

Files ignored due to path filters: 1
  • docs/_generated/rescue/knobs.md is excluded by !**/_generated/**, !docs/_generated/**
Files selected for processing: 22
  • docs/src/developer-guide/rescue-kswv.md
  • docs/src/whats-different/equivalence.md
  • docs/src/whats-different/performance.md
  • src/FMI_search.cpp
  • src/FMI_search.h
  • src/bandedSWA.cpp
  • src/bandedSWA.h
  • src/bandedSWA_compact512.inc
  • src/bntseq.cpp
  • src/bwa.cpp
  • src/bwamem.cpp
  • src/ext_reverse_copy.h
  • src/fast_reader_bseq.c
  • src/fast_reader_bseq.h
  • src/fastmap.cpp
  • src/kswv.cpp
  • src/kswv.h
  • src/main.cpp
  • src/profiling.cpp
  • src/profiling.h
  • src/rescue_env.h
  • src/ungapped_ext.h
  • Autofix · Keep fixing CodeRabbit findings and required CI, and resolving merge conflicts

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands.

@nh13

nh13 commented Oct 28, 2026

Copy link
Copy Markdown

@coderabbitai pause

@coderabbitai

coderabbitai Bot commented Oct 28, 2026

Copy link
Copy Markdown
✅ Action performed

Reviews paused.

@jmschrei

Copy link
Copy Markdown
Author

Hi Nils, hope you're doing well. I felt bad about the way I worded my comment on LinkedIn and so thought I'd see if I could help out here with this nice effort. The above message was written by Claude, but this comment is written by me (though based on my other auto-optimization PRs).

I have been developing a tool, called [auto-optimize], which basically applies the auto-research loop to code to improve speed and memory usage of code, but with extensive robustness checks to try to ensure nothing subtly breaks. I have been applying [auto-optimize] across the genomics ecosystem, starting with my own tools but now extending to other important ecosystem tools, with the hope of making things broadly better. Most recently, I've applied it to MACS3, Bowtie2, chromap, and bedtools.

In an effort to be more robust, the process first augments the unit tests to make sure that nothing breaks during the search. See #654 (feel free to merge that independently of any choice you make here). Any iteration in our loop that broke any unit test was automatically discarded, and we used a set of external validations that were not directly optimized against each iteration to make sure we got identical results. At the end, we ran tens of thousands of cases across assays and file formats, covering every command-line flag, combinations, etc. Importantly, we also created a set of malformed data to make sure that not only did you get identical outputs on good data, but that it failed the same way on bad data so that we didn't just get rid of safeties that were necessary but slow.

If you look at the results below, you'll see that there were many small improvements that, together, led to speed improvements, rather than just a single bottleneck. Your code was already pretty well optimized compared to most things I've run this against, so the gain wasn't too big, but I thought you might be interested in it nonetheless.

image

There are several paths forward here: (1) you can merge the PR as-is, which I'm guessing is unlikely, (2) you can push your own commits on top of this to revert any changes you dislike or build upon them, (3) you can cherry-pick out any individual or set of improvements you like and merge them as separate PRs (as Ben is doing with bowtie2 right now), (4) you can ignore this entirely. We're happy to share any logs or information that would help you!

@nh13

nh13 commented Oct 28, 2026

Copy link
Copy Markdown

Hi @jmschrei, first of all, thank-you for posting here and your thoughtful response, it is most welcome and appreciated.

And also, thank-you for running auto-optimize, I'd welcome a link so I can try it out on a few other tools if that's possible.

For this PR, I am going to leave it open, since it's certainly most welcome. To get your work merged, I am going to split the PR into smaller pieces, and test them across different machine types to see if they also transfer. Then I'll open individual PRs and merge them down, giving you attribution along the way.

Thanks again!

@jmschrei

Copy link
Copy Markdown
Author

Sounds great. Looking forward to seeing the changes included!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants