hmmer-3.3/src/fm_sse.c
/* Function:  fm_getOccCount()
 * Synopsis:  Compute number of occurrences of c in BWT[1..pos]
 *
 * Purpose:   Scan through the BWT to compute number of occurrence of c in BWT[0..pos],
 *            using SSE to scan 16 bytes-at-a-time.
 *
 *            First, use checkpointed occurrence counts in the arrays occCnts_sb and occCnts_b.
 *            The checkpoint used is the one closest to pos, possibly requiring that counts be
 *            subtracted from the checkpointed total
 *
 *            The counting method is SIMD, loading 16 bytes (32 or 64 chars, depending on
 *            alphabet) at a time into the vector co-processors, then counting occurrences. One
 *            constraint of this approach is that occCnts_b checkpoints must be spaced at least
 *            every 32 or 64 chars (16 bytes, in pressed format), and in multiples of 64/32, so
 *            that _mm_load_si128 calls appropriately meet 16-byte-alignment requirements. That's
 *            a reasonable expectation, as spacings of 256 or more seem to give the best speed,
 *            and certainly better space-utilization.
 */
fm_getOccCount()
Synopsis:  Compute number of occurrences of c in BWT[1..pos]

Scan through the BWT to compute number of occurrence of c in BWT[0..pos],
using SSE to scan 16 bytes-at-a-time.

First, use checkpointed occurrence counts in the arrays occCnts_sb and occCnts_b.
The checkpoint used is the one closest to pos, possibly requiring that counts be
subtracted from the checkpointed total

The counting method is SIMD, loading 16 bytes (32 or 64 chars, depending on
alphabet) at a time into the vector co-processors, then counting occurrences. One
constraint of this approach is that occCnts_b checkpoints must be spaced at least
every 32 or 64 chars (16 bytes, in pressed format), and in multiples of 64/32, so
that _mm_load_si128 calls appropriately meet 16-byte-alignment requirements. That's
a reasonable expectation, as spacings of 256 or more seem to give the best speed,
and certainly better space-utilization.
fm_getOccCount () Synopsis : Compute number of occurrences of c in BWT [ 1 .. pos ] Scan through the BWT to compute number of occurrence of c in BWT [ 0 .. pos ] using SSE to scan 16 bytes - at - a - time . First use checkpointed occurrence counts in the arrays occCnts_sb and occCnts_b . The checkpoint used is the one closest to pos possibly requiring that counts be subtracted from the checkpointed total The counting method is SIMD loading 16 bytes ( 32 or 64 chars depending on alphabet ) at a time into the vector co - processors then counting occurrences . One constraint of this approach is that occCnts_b checkpoints must be spaced at least every 32 or 64 chars ( 16 bytes in pressed format ) and in multiples of 64 / 32 so that _mm_load_si128 calls appropriately meet 16 - byte - alignment requirements . That ' s a reasonable expectation as spacings of 256 or more seem to give the best speed and certainly better space - utilization .
int
fm_getOccCount (const FM_DATA *fm, const FM_CFG *cfg, int pos, uint8_t c)
{
  FM_METADATA *meta = cfg->meta;
  int cnt = 0;
  const int b_pos          = (pos+1) / meta->freq_cnt_b ; 
  const int cnt_mod_mask_b = meta->freq_cnt_b - 1; 
  const int b_rel_pos      = (pos+1) & cnt_mod_mask_b; 
  int up_b           = 2*b_rel_pos/meta->freq_cnt_b; 
  int landmark       = ((b_pos+up_b)*meta->freq_cnt_b) - 1 ;
  if (landmark >= fm->N) { 
    up_b      = 0;
    landmark  = (b_pos*(meta->freq_cnt_b)) - 1 ;
  }
#if defined (eslENABLE_SSE)
  int i;
  const uint16_t * occCnts_b  = fm->occCnts_b;
  const uint32_t * occCnts_sb = fm->occCnts_sb;
  const int sb_pos = (pos+1) / meta->freq_cnt_sb; 
  cnt =  FM_OCC_CNT(sb, sb_pos, c );
  if (up_b)
    cnt += FM_OCC_CNT(b, b_pos + 1, c ) ;
  else if ( b_pos !=  sb_pos * (meta->freq_cnt_sb / meta->freq_cnt_b) )
    cnt += FM_OCC_CNT(b, b_pos, c )  ;
  if ( landmark < fm->N || landmark == -1 ) {
    const uint8_t * BWT = fm->BWT;
    register __m128i c_v = *(cfg->fm_chars_v + c);
    register __m128i BWT_v;
    register __m128i tmp_v;
    register __m128i tmp2_v;
    register __m128i counts_v = cfg->fm_neg128_v; 
    if (meta->alph_type == fm_DNA ) {
      if (!up_b) { 
        for (i=1+floor(landmark/4.0) ; i+15<( (pos+1)/4);  i+=16) { 
          BWT_v    = *(__m128i*)(BWT+i);
          FM_MATCH_2BIT(BWT_v, c_v, tmp_v, tmp2_v, tmp_v);
          FM_COUNT_2BIT(tmp_v, tmp2_v, counts_v);
        }
        int remaining_cnt = pos + 1 -  i*4 ;
        if (remaining_cnt > 0) {
          BWT_v    = *(__m128i*)(BWT+i);
          FM_MATCH_2BIT(BWT_v, c_v, tmp_v, tmp2_v, tmp_v);
          tmp_v    = _mm_and_si128(tmp_v, *(cfg->fm_masks_v + remaining_cnt)); 
          FM_COUNT_2BIT(tmp_v, tmp2_v, counts_v);
        }
      } else { 
        for (i=(landmark/4)-15 ; i>(pos/4);  i-=16) {
          BWT_v = *(__m128i*)(BWT+i);
          FM_MATCH_2BIT(BWT_v, c_v, tmp_v, tmp2_v, tmp_v);
          FM_COUNT_2BIT(tmp_v, tmp2_v, counts_v);
        }
        int remaining_cnt = 64 - (pos + 1 - i*4);
        if (remaining_cnt > 0) {
          BWT_v = *(__m128i*)(BWT+i);
          FM_MATCH_2BIT(BWT_v, c_v, tmp_v, tmp2_v, tmp_v);
          tmp_v    = _mm_and_si128(tmp_v, *(cfg->fm_reverse_masks_v + remaining_cnt)); 
          FM_COUNT_2BIT(tmp_v, tmp2_v, counts_v);
        }
      }
/*
    } else if ( meta->alph_type == fm_DNA_full ) {
      if (!up_b) { 
        for (i=1+floor(landmark/2.0) ; i+15<( (pos+1)/2);  i+=16) { 
          BWT_v    = *(__m128i*)(BWT+i);
          FM_MATCH_4BIT(BWT_v, c_v, tmp_v, tmp2_v);
          FM_COUNT_4BIT(tmp_v, tmp2_v, counts_v);
        }
        int remaining_cnt = pos + 1 -  i*2 ;
        if (remaining_cnt > 0) {
          BWT_v    = *(__m128i*)(BWT+i);
          FM_MATCH_4BIT(BWT_v, c_v, tmp_v, tmp2_v);
          tmp_v     = _mm_and_si128(tmp_v, *(cfg->fm_masks_v + (remaining_cnt+1)/2)); 
          tmp2_v    = _mm_and_si128(tmp2_v, *(cfg->fm_masks_v + remaining_cnt/2));
          FM_COUNT_4BIT(tmp_v, tmp2_v, counts_v);
        }
      } else { 
        for (i=(landmark/2)-15 ; i>(pos/2);  i-=16) {
          BWT_v = *(__m128i*)(BWT+i);
          FM_MATCH_4BIT(BWT_v, c_v, tmp_v, tmp2_v);
          FM_COUNT_4BIT(tmp_v, tmp2_v, counts_v);
        }
        int remaining_cnt = 32 - (pos + 1 - i*2);
        if (remaining_cnt > 0) {
          BWT_v = *(__m128i*)(BWT+i);
          FM_MATCH_4BIT(BWT_v, c_v, tmp_v, tmp2_v);
          tmp_v     = _mm_and_si128(tmp_v, *(cfg->fm_reverse_masks_v + remaining_cnt/2)); 
          tmp2_v    = _mm_and_si128(tmp2_v, *(cfg->fm_reverse_masks_v + (remaining_cnt+1)/2));
          FM_COUNT_4BIT(tmp_v, tmp2_v, counts_v);
        }
      }
*/
    } else { 
      if (!up_b) { 
        for (i=1+landmark ; i+15<(pos+1);  i+=16) { 
          BWT_v    = *(__m128i*)(BWT+i);
          BWT_v    = _mm_cmpeq_epi8(BWT_v, c_v);  
          counts_v = _mm_subs_epi8(counts_v, BWT_v); 
        }
        int remaining_cnt = pos + 1 -  i ;
        if (remaining_cnt > 0) {
          BWT_v    = *(__m128i*)(BWT+i);
          BWT_v    = _mm_cmpeq_epi8(BWT_v, c_v);
          BWT_v    = _mm_and_si128(BWT_v, *(cfg->fm_masks_v + remaining_cnt));
          counts_v = _mm_subs_epi8(counts_v, BWT_v);
        }
      } else { 
        for (i=landmark-15 ; i>pos;  i-=16) {
          BWT_v = *(__m128i*)(BWT+i);
          BWT_v    = _mm_cmpeq_epi8(BWT_v, c_v);  
          counts_v = _mm_subs_epi8(counts_v, BWT_v); 
        }
        int remaining_cnt = 16 - (pos + 1 - i);
        if (remaining_cnt > 0) {
          BWT_v = *(__m128i*)(BWT+i);
          BWT_v    = _mm_cmpeq_epi8(BWT_v, c_v);
          BWT_v     = _mm_and_si128(BWT_v, *(cfg->fm_reverse_masks_v + remaining_cnt));
          counts_v = _mm_subs_epi8(counts_v, BWT_v);
        }
      }
    }
    counts_v = _mm_xor_si128(counts_v, cfg->fm_neg128_v); 
    FM_GATHER_8BIT_COUNTS(counts_v,counts_v,counts_v);
    cnt  +=   ( up_b == 1 ?  -1 : 1) * ( _mm_extract_epi16(counts_v, 0) );
  }
  if (c==0 && pos >= fm->term_loc) { 
    cnt--;
  }
#endif 
  return cnt ;
}
int fm_getOccCount ( const FM_DATA * fm , const FM_CFG * cfg , int pos , uint8_t c ) { FM_METADATA * meta = cfg -> meta ; int cnt = 0 ; const int b_pos = ( pos + 1 ) / meta -> freq_cnt_b ; const int cnt_mod_mask_b = meta -> freq_cnt_b - 1 ; const int b_rel_pos = ( pos + 1 ) & cnt_mod_mask_b ; int up_b = 2 * b_rel_pos / meta -> freq_cnt_b ; int landmark = ( ( b_pos + up_b ) * meta -> freq_cnt_b ) - 1 ; if ( landmark >= fm -> N ) { up_b = 0 ; landmark = ( b_pos * ( meta -> freq_cnt_b ) ) - 1 ; } #if defined ( eslENABLE_SSE ) int i ; const uint16_t * occCnts_b = fm -> occCnts_b ; const uint32_t * occCnts_sb = fm -> occCnts_sb ; const int sb_pos = ( pos + 1 ) / meta -> freq_cnt_sb ; cnt = FM_OCC_CNT ( sb , sb_pos , c ) ; if ( up_b ) cnt += FM_OCC_CNT ( b , b_pos + 1 , c ) ; else if ( b_pos != sb_pos * ( meta -> freq_cnt_sb / meta -> freq_cnt_b ) ) cnt += FM_OCC_CNT ( b , b_pos , c ) ; if ( landmark < fm -> N || landmark == -1 ) { const uint8_t * BWT = fm -> BWT ; register __m128i c_v = * ( cfg -> fm_chars_v + c ) ; register __m128i BWT_v ; register __m128i tmp_v ; register __m128i tmp2_v ; register __m128i counts_v = cfg -> fm_neg128_v ; if ( meta -> alph_type == fm_DNA ) { if ( ! up_b ) { for ( i = 1 + floor ( landmark / 4.0 ) ; i + 15 < ( ( pos + 1 ) / 4 ) ; i += 16 ) { BWT_v = * ( __m128i * ) ( BWT + i ) ; FM_MATCH_2BIT ( BWT_v , c_v , tmp_v , tmp2_v , tmp_v ) ; FM_COUNT_2BIT ( tmp_v , tmp2_v , counts_v ) ; } int remaining_cnt = pos + 1 - i * 4 ; if ( remaining_cnt > 0 ) { BWT_v = * ( __m128i * ) ( BWT + i ) ; FM_MATCH_2BIT ( BWT_v , c_v , tmp_v , tmp2_v , tmp_v ) ; tmp_v = _mm_and_si128 ( tmp_v , * ( cfg -> fm_masks_v + remaining_cnt ) ) ; FM_COUNT_2BIT ( tmp_v , tmp2_v , counts_v ) ; } } else { for ( i = ( landmark / 4 ) - 15 ; i > ( pos / 4 ) ; i -= 16 ) { BWT_v = * ( __m128i * ) ( BWT + i ) ; FM_MATCH_2BIT ( BWT_v , c_v , tmp_v , tmp2_v , tmp_v ) ; FM_COUNT_2BIT ( tmp_v , tmp2_v , counts_v ) ; } int remaining_cnt = 64 - ( pos + 1 - i * 4 ) ; if ( remaining_cnt > 0 ) { BWT_v = * ( __m128i * ) ( BWT + i ) ; FM_MATCH_2BIT ( BWT_v , c_v , tmp_v , tmp2_v , tmp_v ) ; tmp_v = _mm_and_si128 ( tmp_v , * ( cfg -> fm_reverse_masks_v + remaining_cnt ) ) ; FM_COUNT_2BIT ( tmp_v , tmp2_v , counts_v ) ; } } } else { if ( ! up_b ) { for ( i = 1 + landmark ; i + 15 < ( pos + 1 ) ; i += 16 ) { BWT_v = * ( __m128i * ) ( BWT + i ) ; BWT_v = _mm_cmpeq_epi8 ( BWT_v , c_v ) ; counts_v = _mm_subs_epi8 ( counts_v , BWT_v ) ; } int remaining_cnt = pos + 1 - i ; if ( remaining_cnt > 0 ) { BWT_v = * ( __m128i * ) ( BWT + i ) ; BWT_v = _mm_cmpeq_epi8 ( BWT_v , c_v ) ; BWT_v = _mm_and_si128 ( BWT_v , * ( cfg -> fm_masks_v + remaining_cnt ) ) ; counts_v = _mm_subs_epi8 ( counts_v , BWT_v ) ; } } else { for ( i = landmark - 15 ; i > pos ; i -= 16 ) { BWT_v = * ( __m128i * ) ( BWT + i ) ; BWT_v = _mm_cmpeq_epi8 ( BWT_v , c_v ) ; counts_v = _mm_subs_epi8 ( counts_v , BWT_v ) ; } int remaining_cnt = 16 - ( pos + 1 - i ) ; if ( remaining_cnt > 0 ) { BWT_v = * ( __m128i * ) ( BWT + i ) ; BWT_v = _mm_cmpeq_epi8 ( BWT_v , c_v ) ; BWT_v = _mm_and_si128 ( BWT_v , * ( cfg -> fm_reverse_masks_v + remaining_cnt ) ) ; counts_v = _mm_subs_epi8 ( counts_v , BWT_v ) ; } } } counts_v = _mm_xor_si128 ( counts_v , cfg -> fm_neg128_v ) ; FM_GATHER_8BIT_COUNTS ( counts_v , counts_v , counts_v ) ; cnt += ( up_b == 1 ? -1 : 1 ) * ( _mm_extract_epi16 ( counts_v , 0 ) ) ; } if ( c == 0 && pos >= fm -> term_loc ) { cnt -- ; } #endif return cnt ; }
Function:  fm_getOccCount()
 Synopsis:  Compute number of occurrences of c in BWT[1..pos]
Function : fm_getOccCount () Synopsis : Compute number of occurrences of c in BWT [ 1 .. pos ]
//floor(pos/b_size)   : the b count element preceding pos //used to compute the mod function // pos % b_size      : how close is pos to the boundary corresponding to b_pos //1 if pos is expected to be closer to the boundary of b_pos+1, 0 otherwise // special case: for a count in the final block, just count from the bottom //floor(pos/sb_size) : the sb count element preceding pos // get the cnt stored at the nearest checkpoint // b_pos has cumulative counts since the prior sb_pos - if sb_pos references the same count as b_pos, it'll doublecount // set to -128, offset to allow each 8-bit int to hold up to 255. // so effectively, can guarantee holding 128*16 = 2048. // Since I count from left or right, whichever is closer, this means // we can support an occ_b interval of up to 4096 with guarantee of // correctness. // count forward, adding // keep running until i begins a run that shouldn't all be counted // leaves only the remaining_cnt chars in the array // count backwards, subtracting // leaves only the remaining_cnt chars in the array /*
    } else if ( meta->alph_type == fm_DNA_full ) {

      if (!up_b) { // count forward, adding
        for (i=1+floor(landmark/2.0) ; i+15<( (pos+1)/2);  i+=16) { // keep running until i begins a run that shouldn't all be counted
          BWT_v    = *(__m128i*)(BWT+i);
          FM_MATCH_4BIT(BWT_v, c_v, tmp_v, tmp2_v);
          FM_COUNT_4BIT(tmp_v, tmp2_v, counts_v);
        }
        int remaining_cnt = pos + 1 -  i*2 ;
        if (remaining_cnt > 0) {
          BWT_v    = *(__m128i*)(BWT+i);
          FM_MATCH_4BIT(BWT_v, c_v, tmp_v, tmp2_v);
          tmp_v     = _mm_and_si128(tmp_v, *(cfg->fm_masks_v + (remaining_cnt+1)/2)); // mask characters we don't want to count
          tmp2_v    = _mm_and_si128(tmp2_v, *(cfg->fm_masks_v + remaining_cnt/2));
          FM_COUNT_4BIT(tmp_v, tmp2_v, counts_v);
        }

      } else { // count backwards, subtracting
        for (i=(landmark/2)-15 ; i>(pos/2);  i-=16) {
          BWT_v = *(__m128i*)(BWT+i);
          FM_MATCH_4BIT(BWT_v, c_v, tmp_v, tmp2_v);
          FM_COUNT_4BIT(tmp_v, tmp2_v, counts_v);
        }

        int remaining_cnt = 32 - (pos + 1 - i*2);
        if (remaining_cnt > 0) {
          BWT_v = *(__m128i*)(BWT+i);
          FM_MATCH_4BIT(BWT_v, c_v, tmp_v, tmp2_v);
          tmp_v     = _mm_and_si128(tmp_v, *(cfg->fm_reverse_masks_v + remaining_cnt/2)); // mask characters we don't want to count
          tmp2_v    = _mm_and_si128(tmp2_v, *(cfg->fm_reverse_masks_v + (remaining_cnt+1)/2));
          FM_COUNT_4BIT(tmp_v, tmp2_v, counts_v);
        }
      }
*/ //amino // count forward, adding // keep running until i begins a run that shouldn't all be counted // each byte is all 1s if matching, all zeros otherwise // adds 1 for each matching byte  (subtracting negative 1) // mask characters we don't want to count // count backwards, subtracting // each byte is all 1s if matching, all zeros otherwise // adds 1 for each matching byte  (subtracting negative 1) // mask characters we don't want to count //tmp2_v    = _mm_and_si128(tmp2_v, *(cfg->fm_reverse_masks_v + (remaining_cnt+1)/2)); //counts are stored in signed bytes, base -128. Move them to unsigned bytes // I overcounted 'A' by one, because '$' was replaced with an 'A' //#if   defined (eslENABLE_SSE)