Refactor SIMD code and add ARM Neon nibble2base() implementation by jmarshall · Pull Request #1795 · samtools/htslib

jmarshall · 2024-06-24T13:28:07Z

This PR moves nibble2base_ssse3() from an internal header file to a new source file, simd.c, and adds an ARM implementation nibble2base_neon() alongside. Because these functions are always called via a pointer, they are never inlined so there is nothing gained (other than build complexity!) by having them in a header.

I initially translated the SSSE3 algorithm directly to Neon, during which I noticed some tweaks that could be usefully applied to the SSSE3 version — see the “Minor improvements” commit below. I then realised that the Neon implementation could be simplified by using Neon's rather excellent interleaved load/store instructions, which I then replaced by an in-register zip as it is (sadly) faster (on aarch64 at least). (It may be possible to apply a similar approach for Intel via the unpack instruction.)

Indicative timings for various versions of the Neon implementation, on ARM and, for interest, an iteration count on x86_64 that makes the nibble2base_default time comparable — numbers are comparable vertically, but not directly horizontally:

Apple M2  i5 SSSE3
 seconds   seconds

 262.389   916.215   nibble2base_single (i.e., scalar code doing one nucleotide at a time)
 186.283   186.778   nibble2base_default (i.e., scalar code doing two nucleotides at a time)
  36.053    46.700   Ruben's nibble2base_ssse3 algorithm
  24.783             interleaved using vst2q_u8
  19.287             interleaved using vzipq_u8

I moved the existing code to the new source file in a pair of commits marked as authored by @rhpvorderman and @daviesrob, so that git blame preserves the authorship of their untouched lines. Hence ideally this PR would be merged via a merge commit rather than squashed.

This function is always called via a function pointer, hence it is never inlined so need not be in a header. [Commit marked as authored by RV so git blame attributes the moved code to Ruben, as it is unchanged from the original code. Temporarily include simd.c from the header so this commit builds and runs. -JM]

Also add copyright boilerplate to the new file. [Commit marked as authored by RMD so git blame attributes the boilerplate and moved code to Rob, as the latter is unchanged from the original code. -JM]

… file Build simd.c as a source file, so there is one copy of nibble2base_ssse3() rather than potentially one per translation unit that uses it. When everything was static, the names needed no prefix. However now that there is a non-static function pointer, it needs an htslib-specific prefix to avoid polluting the namespace. (In libhts.so it will usually be a non-visible symbol, but in libhts.a it needs a prefix.)

Configure checks for __builtin_cpu_supports("ssse3"), so the check for __attribute__((target)) might as well be there too. This enables this internal-only infrastructure to be removed from htslib/hts_defs.h, which is primarily for defines needed in the public headers. Add a configure check for __attribute__((constructor)) too, as we will likely be using it on other SIMD-using platforms as well.

Notably the adjustment to seq_vec_end_ptr stops it from reverting to the scalar code prematurely: at present it drops back to scalar when <=32 entries remain, rather than <=31.

As ARM CPU capability instructions are privileged, we also need to add cpu_supports_neon() and implement it to query for Neon/AdvSIMD in various platform-dependent ways. (On 32-bit ARM, glibc helpfully renamed the HWCAP_* macros to HWCAP_ARM_*, so accept both on Linux.) 32-bit ARM GCC erroneously does not define vst1q_u8_x2() (bug 71233, fixed in v14.1 in this case), so avoid it on 32-bit ARM by writing interleaved via vst2q_u8() instead.

jkbonfield · 2024-06-24T14:24:22Z

That's impressive. On my intel machine SSSE3 appears to be running about 4x faster than the default version. You're getting 8-9x speed up which is great.

PS. The mythical htslib v2.0 will not use nibble encoding in memory. It's a royal pain in the derriere!

jmarshall · 2024-06-25T11:54:39Z

I updated the timings to values from the final code and test code in this PR (previously the timings were from earlier development code), included the 1-nibble-at-a-time naive scalar code timings, and added some x86_64 timings on a different iteration count that makes the times on one row comparable. (Not totally sure I believe that 916 seconds! And I only ran it once. But it's indicative that the difference between the two scalar versions is much bigger on x86_64.)

It's not really meaningful to compare between the columns, but what it does tell you is this: the 2-nibbles-at-a-time scalar code is a much smaller win over 1-nibble-at-a-time on ARM than on x86_64. So this suggests that the apparent 8-9x speed up might not be that astonishing, because there is a lot left on the table by the 2-at-a-time scalar code on ARM, at least as it is written here. Or something — who knows, it's a nice speed up anyway 😄

jkbonfield · 2024-06-25T13:50:47Z

Yep, very nice and thank you.

I hunted the origin of the nibble2base and backtracked it to my initial BAM implementation in io_lib in Jan 2013 :-)
It's since been tarted up a bit and avoiding the need for unaligned access ifdefs by use of memcpy, but essentially the same deal. It got retrofitted into the BAM handling code in htslib shortly after CRAM was added.

https://github.com/jkbonfield/io_lib/blob/master/io_lib/bam.c#L3650-L3657

rhpvorderman · 2024-06-29T19:47:56Z

It may be possible to apply a similar approach for Intel via the unpack instruction.

I think it would be possible, and possibly more optimal. Interesting find!

rhpvorderman · 2024-06-29T20:03:17Z

Did a quick test in Sequali, the following code is simpler and all tests pass:

static void 
decode_bam_sequence_ssse3(uint8_t *dest, const uint8_t *encoded_sequence, size_t length) 
{

    static const uint8_t *nuc_lookup = (uint8_t *)"=ACMGRSVTWYHKDBN";
    const uint8_t *dest_end_ptr = dest + length;
    uint8_t *dest_cursor = dest;
    const uint8_t *encoded_cursor = encoded_sequence;
    const uint8_t *dest_vec_end_ptr = dest_end_ptr - (2 * sizeof(__m128i) - 1);
    __m128i nuc_lookup_vec = _mm_lddqu_si128((__m128i *)nuc_lookup);
    /* Nucleotides are encoded 4-bits per nucleotide and stored in 8-bit bytes 
       as follows: |AB|CD|EF|GH|. The 4-bit codes (going from 0-15) can be used 
       together with the pshufb instruction as a lookup table. The most efficient
       way is to use bitwise AND and shift to create two vectors. One with all 
       the upper codes (|A|C|E|G|) and one with the lower codes (|B|D|F|H|). 
       The lookup can then be performed and the resulting vectors can be 
       interleaved again using the unpack instructions. */
    while (dest_cursor < dest_vec_end_ptr) {
        __m128i encoded = _mm_lddqu_si128((__m128i *)encoded_cursor);
        __m128i encoded_upper = _mm_srli_epi64(encoded, 4);
        encoded_upper = _mm_and_si128(encoded_upper, _mm_set1_epi8(15));
        __m128i encoded_lower = _mm_and_si128(encoded, _mm_set1_epi8(15));
        __m128i nucs_upper = _mm_shuffle_epi8(nuc_lookup_vec, encoded_upper);
        __m128i nucs_lower = _mm_shuffle_epi8(nuc_lookup_vec, encoded_lower);
        __m128i first_nucleotides = _mm_unpacklo_epi8(nucs_upper, nucs_lower);
        __m128i second_nucleotides = _mm_unpackhi_epi8(nucs_upper, nucs_lower);
        _mm_storeu_si128((__m128i *)dest_cursor, first_nucleotides);
        _mm_storeu_si128((__m128i *)(dest_cursor + 16), second_nucleotides);
        encoded_cursor += sizeof(__m128i);
        dest_cursor += 2 * sizeof(__m128i);
    }
    decode_bam_sequence_default(dest_cursor, encoded_cursor, dest_end_ptr - dest_cursor);
}

Thanks @jmarshall! This is a great improvement codewise!

EDIT: (Of course I will contribute this after this is merged, or feel free to incorporate it right away).

daviesrob · 2024-07-01T11:21:31Z

When building on my Mac without using ./configure first, I get:

simd.c:63:19: warning: unused function 'cpu_supports_neon' [-Wunused-function]
static inline int cpu_supports_neon(void) {

This is because use of the function is gated by defined __ARM_NEON but its definition is not. Presumably this warning would also appear when building on an ARM variant that doesn't have NEON.

jmarshall · 2024-07-01T12:02:31Z

Oops. I think it was more because on ARM non-configure builds, it wasn't building nibble2base_neon() because HAVE_ATTRIBUTE_CONSTRUCTOR wasn't defined. Fixed now.

Probably in theory the definition of cpu_supports_neon() should be gated by #if defined BUILDING_SIMD_NIBBLE2BASE || defined BUILDING_SIMD_WHATEVER_ELSE || … but life is short and ARM compilers should be supporting building these functions even if the eventual worker CPU doesn't support running them…

daviesrob · 2024-07-01T13:51:00Z

Cheers, that's fixed the problem.

Notice: this is the last SAMtools / HTSlib release where CRAM 3.0 will be the default CRAM version. From the next we will change to CRAM 3.1 unless the version is explicitly specified, for example using "samtools view -O cram,version=3.0". Updates ------- * Extend annot-tsv with several new command line options. --delim permits use of other delimiters. --headers for selection of other header formats. --no-header-idx to suppress column index numbers in header. Also removed -h as it is now short for --headers. Note --help still works. (PR samtools#1779) * Allow annot-tsv -a to rename annotations. (PR samtools#1709) * Extend annot-tsv --overlap to be able to specify the overlap fraction separately for source and target. (PR samtools#1811) * Added new APIs to facilitate low-level CRAM container manipulations, used by the new "samtools cat" region filtering code. Functions are: cram_container_get_coords() cram_filter_container() cram_index_extents() cram_container_num2offset() cram_container_offset2num() cram_num_containers() cram_num_containers_between() Also improved cram_index_query() to cope with HTS_IDX_NOCOOR regions. (PR samtools#1771) * Bgzip now retains file modification and access times when compressing and decompressing. (PR samtools#1727, fixes samtools#1718. Requested by Gert Hulselmans.) * Use FNV1a for string hashing in khash. The old algorithm was particularly weak with base-64 style strings and lead to a large number of collisions. (PR samtools#1806. Fixes samtools/samtools#2066, reported by Hans-Joachim Ruscheweyh) * Improve the speed of the nibble2base() function on Intel (PR samtools#1667, PR samtools#1764, PR samtools#1786, PR samtools#1802, thanks to Ruben Vorderman) and ARM (PR samtools#1795, thanks to John Marshall). * bgzf_getline() will now warn if it encounters UTF-16 data. (PR samtools#1487, thanks to John Marshall) * Speed up bgzf_read(). While this does not reduce CPU significantly, it does increase the maximum parallelism available permitting 10-15% faster decoding. (PR samtools#1772, PR samtools#1800, Issue samtools#1798) * Speed up faidx by use of better isgraph methods (PR samtools#1797) and whole-line reading (PR samtools#1799, thanks to John Marshall). * Speed up kputll() function, speeding up BAM -> SAM conversion by about 5% and also samtools depth. (PR samtools#1805) * Added more example code, covering fasta/fastq indexing, tabix indexing and use of the thread pool. (PR samtools#1666) Build Changes ------------- * Code warning fixes for pedantic compilers (PR samtools#1777) and avoid some undefined behaviour (PR samtools#1810, PR samtools#1816, PR samtools#1828). * Windows based CI has been migrated from AppVeyor to GitHub Actions. (PR samtools#1796, PR samtools#1803, PR samtools#1808) * Miscellaneous minor build infrastructure and code fixes. (PR samtools#1807, PR samtools#1829, both thanks to John Marshall) * Updated htscodecs submodule to version 1.6.1 (PR samtools#1828) * Fixed an awk script in the Makefile that only worked with gawk. (PR samtools#1831) Bug fixes --------- * Fix small OSS-Fuzz reported issues with CRAM encoding and long CIGARS and/or illegal positions. (PR samtools#1775, PR samtools#1801, PR samtools#1817) * Fix issues with on-the-fly indexing of VCF/BCF (bcftools --write-index) when not using multiple threads. (PR samtools#1837. Fixes samtools/bcftools#2267, reported by Giulio Genovese) * Stricter limits on POS / MPOS / TLEN in sam_parse1(). This fixes a signed overflow reported by OSS-Fuzz and should help prevent other as-yet undetected bugs. (PR samtools#1812) * Check that the underlying file open worked for preload: URLs. Fixes a NULL pointer dereference reported by OSS-Fuzz. (PR samtools#1821) * Fix an infinite loop in hts_itr_query() when given extremely large positions which cause integer overflow. Also adds hts_bin_maxpos() and hts_idx_maxpos() functions. (PR samtools#1774, thanks to John Marshall and reported by Jesus Alberto Munoz Mesa) * Fix an out of bounds read in hts_itr_multi_next() when switching chromosomes. This bug is present in releases 1.11 to 1.20. (PR samtools#1788. Fixes samtools/samtools#2063, reported by acorvelo) * Work around parsing problems with colons in CHROM names. Fixes samtools/bcftools#2139. (PR samtools#1781, John Marshall / James Bonfield) * Correct the CPU detection for Mac OS X 10.7. cpuid is used by htscodecs (see samtools/htscodecs#116), and the corresponding changes in htslib are PR samtools#1785. Reported by Ryan Carsten Schmidt. * Make BAM zero-length intervals work the same as CRAM; permitted and returning overlapping records. (PR samtools#1787. Fixes samtools/samtools#2060, reported by acorvelo) * Replace assert() with abort() in BCF synced reader. This is not an ideal solution, but it gives consistent behaviour when compiling with or without NDEBUG. (PR samtools#1791, thanks to Martin Pollard) * Fixed failure to change the write block size on compressed SAM or VCF files due to an internal type confusion. (PR samtools#1826) * Fixed an out-of-bounds read in cram_codec_iter_next() (PR samtools#1832)

htslib release 1.21: The primary user-visible changes in this release are updates to the annot-tsv tool and some speed improvements. Full details of other changes and bugs fixed are below. Notice: this is the last SAMtools / HTSlib release where CRAM 3.0 will be the default CRAM version. From the next we will change to CRAM 3.1 unless the version is explicitly specified, for example using "samtools view -O cram,version=3.0". Updates ------- * Extend annot-tsv with several new command line options. --delim permits use of other delimiters. --headers for selection of other header formats. --no-header-idx to suppress column index numbers in header. Also removed -h as it is now short for --headers. Note --help still works. (PR samtools#1779) * Allow annot-tsv -a to rename annotations. (PR samtools#1709) * Extend annot-tsv --overlap to be able to specify the overlap fraction separately for source and target. (PR samtools#1811) * Added new APIs to facilitate low-level CRAM container manipulations, used by the new "samtools cat" region filtering code. Functions are: cram_container_get_coords() cram_filter_container() cram_index_extents() cram_container_num2offset() cram_container_offset2num() cram_num_containers() cram_num_containers_between() Also improved cram_index_query() to cope with HTS_IDX_NOCOOR regions. (PR samtools#1771) * Bgzip now retains file modification and access times when compressing and decompressing. (PR samtools#1727, fixes samtools#1718. Requested by Gert Hulselmans.) * Use FNV1a for string hashing in khash. The old algorithm was particularly weak with base-64 style strings and lead to a large number of collisions. (PR samtools#1806. Fixes samtools/samtools#2066, reported by Hans-Joachim Ruscheweyh) * Improve the speed of the nibble2base() function on Intel (PR samtools#1667, PR samtools#1764, PR samtools#1786, PR samtools#1802, thanks to Ruben Vorderman) and ARM (PR samtools#1795, thanks to John Marshall). * bgzf_getline() will now warn if it encounters UTF-16 data. (PR samtools#1487, thanks to John Marshall) * Speed up bgzf_read(). While this does not reduce CPU significantly, it does increase the maximum parallelism available permitting 10-15% faster decoding. (PR samtools#1772, PR samtools#1800, Issue samtools#1798) * Speed up faidx by use of better isgraph methods (PR samtools#1797) and whole-line reading (PR samtools#1799, thanks to John Marshall). * Speed up kputll() function, speeding up BAM -> SAM conversion by about 5% and also samtools depth. (PR samtools#1805) * Added more example code, covering fasta/fastq indexing, tabix indexing and use of the thread pool. (PR samtools#1666) Build Changes ------------- * Code warning fixes for pedantic compilers (PR samtools#1777) and avoid some undefined behaviour (PR samtools#1810, PR samtools#1816, PR samtools#1828). * Windows based CI has been migrated from AppVeyor to GitHub Actions. (PR samtools#1796, PR samtools#1803, PR samtools#1808) * Miscellaneous minor build infrastructure and code fixes. (PR samtools#1807, PR samtools#1829, both thanks to John Marshall) * Updated htscodecs submodule to version 1.6.1 (PR samtools#1828) * Fixed an awk script in the Makefile that only worked with gawk. (PR samtools#1831) Bug fixes --------- * Fix small OSS-Fuzz reported issues with CRAM encoding and long CIGARS and/or illegal positions. (PR samtools#1775, PR samtools#1801, PR samtools#1817) * Fix issues with on-the-fly indexing of VCF/BCF (bcftools --write-index) when not using multiple threads. (PR samtools#1837. Fixes samtools/bcftools#2267, reported by Giulio Genovese) * Stricter limits on POS / MPOS / TLEN in sam_parse1(). This fixes a signed overflow reported by OSS-Fuzz and should help prevent other as-yet undetected bugs. (PR samtools#1812) * Check that the underlying file open worked for preload: URLs. Fixes a NULL pointer dereference reported by OSS-Fuzz. (PR samtools#1821) * Fix an infinite loop in hts_itr_query() when given extremely large positions which cause integer overflow. Also adds hts_bin_maxpos() and hts_idx_maxpos() functions. (PR samtools#1774, thanks to John Marshall and reported by Jesus Alberto Munoz Mesa) * Fix an out of bounds read in hts_itr_multi_next() when switching chromosomes. This bug is present in releases 1.11 to 1.20. (PR samtools#1788. Fixes samtools/samtools#2063, reported by acorvelo) * Work around parsing problems with colons in CHROM names. Fixes samtools/bcftools#2139. (PR samtools#1781, John Marshall / James Bonfield) * Correct the CPU detection for Mac OS X 10.7. cpuid is used by htscodecs (see samtools/htscodecs#116), and the corresponding changes in htslib are PR samtools#1785. Reported by Ryan Carsten Schmidt. * Make BAM zero-length intervals work the same as CRAM; permitted and returning overlapping records. (PR samtools#1787. Fixes samtools/samtools#2060, reported by acorvelo) * Replace assert() with abort() in BCF synced reader. This is not an ideal solution, but it gives consistent behaviour when compiling with or without NDEBUG. (PR samtools#1791, thanks to Martin Pollard) * Fixed failure to change the write block size on compressed SAM or VCF files due to an internal type confusion. (PR samtools#1826) * Fixed an out-of-bounds read in cram_codec_iter_next() (PR samtools#1832) # -----BEGIN PGP SIGNATURE----- # # iQJJBAABCgAzFiEEaGn1wQ1nqxSLDC2XHsvluhXGTMwFAmbjBlMVHHJtZCtnaXRA # c2FuZ2VyLmFjLnVrAAoJEB7L5boVxkzMgb0P/3AAxkJNLGK8Q1pBBo+in7KzwBOc # /nRfvReki47r4+i1ev98abEF/+XZoYSS8ZkiJH8JMkLCclb1ch2fDQTIL/sjJwrU # MfhYx6oHhzLUOjpG24m390wlvlkEsg5S69jbgImWxLQ1n4ZcSuZGumiwPx5274w5 # 3R/K9vIXtuSF9DTD/jrTgoZeBiN+B4XefAJqHLSo4/XzDo7DNMSZryNGFpHdyKcB # oEfEXX45vJLgzA3Fxo3LqNpC8UTDO5C9u9Uup+hCLN7RTBjV5Fu0UxcwWTSPqfiC # Jroi2FhxNPlDkRL3jraI8LrB3Cqb0G/31OU+cryasrqYnfJrl0e/MyXiRSa0sg25 # uHHWVyyC4JR7uBfUCppTVf76iajbkIVpHgHLZS1heGwM4PCL+AalvsIgGcvLXcBs # g/AwZLoXzao+pbzg8hzN1JuLB5hHcrHp6TsQ6JbTTrWaLeAL3s/+eQ4+d+Np2iUn # v8TMg6onu1HiNgVMQcFXSLC2Q51rB+lD6AA3eusKjcQsaxyk+U7SLCx3OhbvxGsd # SjxQWP1W+WW//e9nm5f4CTSSf1PHoDMyR/kVEc3qwEWJslri68RyB69HrDIVzrWT # 3gKAKJYafi5ssFHsKTsCngESpDsTg3aQV1AGOQPLN/TIjyYDTUOL+/HylRs3k9XN # KOGO5yggx+yqeL3w # =lHVl # -----END PGP SIGNATURE----- # gpg: Signature made Thu 12 Sep 2024 08:18:43 AM PDT # gpg: using RSA key 6869F5C10D67AB148B0C2D971ECBE5BA15C64CCC # gpg: issuer "rmd+git@sanger.ac.uk" # gpg: Can't check signature: No public key

rhpvorderman and others added 7 commits June 2, 2024 19:32

Move nibble2base_resolve() infrastructure to simd.c

87186f8

Also add copyright boilerplate to the new file. [Commit marked as authored by RMD so git blame attributes the boilerplate and moved code to Rob, as the latter is unchanged from the original code. -JM]

Minor improvements to Intel implementation

a5ba1a6

Notably the adjustment to seq_vec_end_ptr stops it from reverting to the scalar code prematurely: at present it drops back to scalar when <=32 entries remain, rather than <=31.

Add correctness and speed tests for nibble2base() routines

2c48f9f

daviesrob assigned jkbonfield Jun 27, 2024

This was referenced Jun 30, 2024

Improve BAM nibble decoding and dispatch functions at load time. rhpvorderman/sequali#177

Merged

Decode BAM sequence with less code. Dispatch function at dynamic link time rather than at run time. marcelm/dnaio#135

Merged

Define HAVE_ATTRIBUTE_CONSTRUCTOR on ARM too for non-configure builds

552905d

jkbonfield merged commit 1b4cda6 into samtools:develop Jul 4, 2024

jmarshall deleted the simd branch July 4, 2024 09:44

rhpvorderman mentioned this pull request Jul 8, 2024

Simplify SSSE3 nibble2base function #1802

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Refactor SIMD code and add ARM Neon nibble2base() implementation#1795

Refactor SIMD code and add ARM Neon nibble2base() implementation#1795
jkbonfield merged 8 commits into
samtools:developfrom
jmarshall:simd

jmarshall commented Jun 24, 2024 •

edited

Loading

Uh oh!

jkbonfield commented Jun 24, 2024

Uh oh!

jmarshall commented Jun 25, 2024 •

edited

Loading

Uh oh!

jkbonfield commented Jun 25, 2024

Uh oh!

rhpvorderman commented Jun 29, 2024

Uh oh!

rhpvorderman commented Jun 29, 2024 •

edited

Loading

Uh oh!

daviesrob commented Jul 1, 2024

Uh oh!

jmarshall commented Jul 1, 2024 •

edited

Loading

Uh oh!

daviesrob commented Jul 1, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Uh oh!

Conversation

jmarshall commented Jun 24, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jkbonfield commented Jun 24, 2024

Uh oh!

jmarshall commented Jun 25, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jkbonfield commented Jun 25, 2024

Uh oh!

rhpvorderman commented Jun 29, 2024

Uh oh!

rhpvorderman commented Jun 29, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

daviesrob commented Jul 1, 2024

Uh oh!

jmarshall commented Jul 1, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

daviesrob commented Jul 1, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

jmarshall commented Jun 24, 2024 •

edited

Loading

jmarshall commented Jun 25, 2024 •

edited

Loading

rhpvorderman commented Jun 29, 2024 •

edited

Loading

jmarshall commented Jul 1, 2024 •

edited

Loading