MikkelSchubert / adapterremoval / #39

Committed 20 Aug 2024 01:31PM UTC coverage: 79.602% (-3.8%) from 83.365%

Build # #39

Build Type

push

travis-ci

Committed by

MikkelSchubert

Commit Message

update schema URL to use github pages

Run Details

2279 of 2863 relevant lines covered (79.6%)

14257.45 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/simd_avx512bw.cpp

/*************************************************************************\
 * AdapterRemoval - cleaning next-generation sequencing reads            *
 *                                                                       *
 * Copyright (C) 2024 by Mikkel Schubert - mikkelsch@gmail.com           *
 *                                                                       *
 * This program is free software: you can redistribute it and/or modify  *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * This program is distributed in the hope that it will be useful,       *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have received a copy of the GNU General Public License     *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>. *
\*************************************************************************/
#include "simd.hpp"    // declarations
#include <bitset>      // for bitset
#include <cstddef>     // for size_t
#include <immintrin.h> // for _mm512_set1_epi8, __m512i, _mm512_lo...

#ifdef AR_SUPPORTS_AVX512

namespace adapterremoval {

namespace simd {

namespace {

/** Counts the number of masked bytes **/
inline auto
count_masked_avx512(__mmask64 value)
{
  // Generate 64 bit mask from most significant bits of each byte and count bits
  return std::bitset<64>(value).count();
}

} // namespace

bool
compare_subsequences_avx512(size_t& n_mismatches,
                            size_t& n_ambiguous,
                            const char* seq_1,
                            const char* seq_2,
                            size_t length,
                            size_t max_penalty)
{
  //! Mask of all Ns
  const auto n_mask = _mm512_set1_epi8('N');

  while (length) {
    __m512i s1;
    __m512i s2;

    if (length < 64) {
      const auto mask = 0xFFFFFFFFFFFFFFFFLLU >> (64 - length);
      s1 = _mm512_maskz_loadu_epi8(mask, seq_1);
      s2 = _mm512_maskz_loadu_epi8(mask, seq_2);

      length = 0;
    } else {
      s1 = _mm512_loadu_epi8(seq_1);
      s2 = _mm512_loadu_epi8(seq_2);

      seq_1 += 64;
      seq_2 += 64;
      length -= 64;
    }

    // Sets 0xFF for every byte where one or both nts is N
    const auto ns_mask =
      _mm512_cmpeq_epu8_mask(s1, n_mask) | _mm512_cmpeq_epu8_mask(s2, n_mask);

    // Sets 0xFF for every byte where bytes are equal or N
    const auto eq_mask = _mm512_cmpeq_epu8_mask(s1, s2) | ns_mask;

    n_mismatches += 64 - count_masked_avx512(eq_mask);
    if (2 * n_mismatches + n_ambiguous > max_penalty) {
      return false;
    }

    // Early termination is almost always due to mismatches, so updating the
    // number of Ns after the above check saves time in the common case.
    n_ambiguous += count_masked_avx512(ns_mask);
  }

  return true;
}

} // namespace simd

} // namespace adapterremoval

#endif

1	/*************************************************************************\
2	* AdapterRemoval - cleaning next-generation sequencing reads *
3	* *
4	* Copyright (C) 2024 by Mikkel Schubert - mikkelsch@gmail.com *
5	* *
6	* This program is free software: you can redistribute it and/or modify *
7	* it under the terms of the GNU General Public License as published by *
8	* the Free Software Foundation, either version 3 of the License, or *
9	* (at your option) any later version. *
10	* *
11	* This program is distributed in the hope that it will be useful, *
12	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
13	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14	* GNU General Public License for more details. *
15	* *
16	* You should have received a copy of the GNU General Public License *
17	* along with this program. If not, see <http://www.gnu.org/licenses/>. *
18	\*************************************************************************/
19	#include "simd.hpp" // declarations
20	#include <bitset> // for bitset
21	#include <cstddef> // for size_t
22	#include <immintrin.h> // for _mm512_set1_epi8, __m512i, _mm512_lo...
23
24	#ifdef AR_SUPPORTS_AVX512
25
26	namespace adapterremoval {
27
28	namespace simd {
29
30	namespace {
31
32	/ Counts the number of masked bytes /
33	inline auto
34	count_masked_avx512(__mmask64 value)	×
35	{
36	// Generate 64 bit mask from most significant bits of each byte and count bits
37	return std::bitset<64>(value).count();	×
38	}
39
40	} // namespace
41
42	bool
43	compare_subsequences_avx512(size_t& n_mismatches,	×
44	size_t& n_ambiguous,
45	const char* seq_1,
46	const char* seq_2,
47	size_t length,
48	size_t max_penalty)
49	{
50	//! Mask of all Ns
51	const auto n_mask = _mm512_set1_epi8('N');	×
52
53	while (length) {	×
54	__m512i s1;	×
55	__m512i s2;	×
56
57	if (length < 64) {	×
58	const auto mask = 0xFFFFFFFFFFFFFFFFLLU >> (64 - length);	×
59	s1 = _mm512_maskz_loadu_epi8(mask, seq_1);	×
60	s2 = _mm512_maskz_loadu_epi8(mask, seq_2);	×
61
62	length = 0;	×
63	} else {
64	s1 = _mm512_loadu_epi8(seq_1);	×
65	s2 = _mm512_loadu_epi8(seq_2);	×
66
67	seq_1 += 64;	×
68	seq_2 += 64;	×
69	length -= 64;	×
70	}
71
72	// Sets 0xFF for every byte where one or both nts is N
73	const auto ns_mask =	×
74	_mm512_cmpeq_epu8_mask(s1, n_mask) \| _mm512_cmpeq_epu8_mask(s2, n_mask);	×
75
76	// Sets 0xFF for every byte where bytes are equal or N
77	const auto eq_mask = _mm512_cmpeq_epu8_mask(s1, s2) \| ns_mask;	×
78
79	n_mismatches += 64 - count_masked_avx512(eq_mask);	×
80	if (2 * n_mismatches + n_ambiguous > max_penalty) {	×
81	return false;
82	}
83
84	// Early termination is almost always due to mismatches, so updating the
85	// number of Ns after the above check saves time in the common case.
86	n_ambiguous += count_masked_avx512(ns_mask);	×
87	}
88
89	return true;
90	}
91
92	} // namespace simd
93
94	} // namespace adapterremoval
95
96	#endif

MikkelSchubert / adapterremoval / #39

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous