text: /**
   * @file
   *
   * @author jeff.daily@pnnl.gov
   *
   * Copyright (c) 2015 Battelle Memorial Institute.
   */
  #include "config.h"
  
  #include <stdint.h>
  #include <stdlib.h>
  #include <string.h>
  
  #include <emmintrin.h>
  #include <smmintrin.h>
  
  #include "parasail.h"
  #include "parasail/memory.h"
  #include "parasail/internal_sse.h"
  
  #define FASTSTATS
  
  #define NEG_INF (INT64_MIN/(int64_t)(2))
  
  static inline __m128i _mm_cmpgt_epi64_rpl(__m128i a, __m128i b) {
      __m128i_64_t A;
      __m128i_64_t B;
      A.m = a;
      B.m = b;
      A.v[0] = (A.v[0]>B.v[0]) ? 0xFFFFFFFFFFFFFFFF : 0;
      A.v[1] = (A.v[1]>B.v[1]) ? 0xFFFFFFFFFFFFFFFF : 0;
      return A.m;
  }
  
  static inline __m128i _mm_max_epi64_rpl(__m128i a, __m128i b) {
      __m128i_64_t A;
      __m128i_64_t B;
      A.m = a;
      B.m = b;
      A.v[0] = (A.v[0]>B.v[0]) ? A.v[0] : B.v[0];
      A.v[1] = (A.v[1]>B.v[1]) ? A.v[1] : B.v[1];
      return A.m;
  }
  
  static inline __m128i _mm_cmplt_epi64_rpl(__m128i a, __m128i b) {
      __m128i_64_t A;
      __m128i_64_t B;
      A.m = a;
      B.m = b;
      A.v[0] = (A.v[0]<B.v[0]) ? 0xFFFFFFFFFFFFFFFF : 0;
      A.v[1] = (A.v[1]<B.v[1]) ? 0xFFFFFFFFFFFFFFFF : 0;
      return A.m;
  }
  
  static inline int64_t _mm_hmax_epi64_rpl(__m128i a) {
      a = _mm_max_epi64_rpl(a, _mm_srli_si128(a, 8));
      return _mm_extract_epi64(a, 0);
  }
  
  
  #ifdef PARASAIL_TABLE
  static inline void arr_store_si128(
          int *array,
          __m128i vH,
          int32_t t,
          int32_t seglen,
          int32_t d,
          int32_t dlen)
  {
      array[(0*seglen+t)*dlen + d] = (int64_t)_mm_extract_epi64(vH, 0);
      array[(1*seglen+t)*dlen + d] = (int64_t)_mm_extract_epi64(vH, 1);
  }
  #endif
  
  #ifdef PARASAIL_ROWCOL
  static inline void arr_store_col(
          int *col,
          __m128i vH,
          int32_t t,
          int32_t seglen)
  {
      col[0*seglen+t] = (int64_t)_mm_extract_epi64(vH, 0);
      col[1*seglen+t] = (int64_t)_mm_extract_epi64(vH, 1);
  }
  #endif
  
  #ifdef PARASAIL_TABLE
  #define FNAME parasail_sw_stats_table_striped_sse41_128_64
  #define PNAME parasail_sw_stats_table_striped_profile_sse41_128_64
  #define INAME PNAME
  #define STATIC
  #else
  #ifdef PARASAIL_ROWCOL
  #define FNAME parasail_sw_stats_rowcol_striped_sse41_128_64
  #define PNAME parasail_sw_stats_rowcol_striped_profile_sse41_128_64
  #define INAME PNAME
  #define STATIC
  #else
  #define FNAME parasail_sw_stats_striped_sse41_128_64
  #ifdef FASTSTATS
  #define PNAME parasail_sw_stats_striped_profile_sse41_128_64_internal
  #define INAME parasail_sw_stats_striped_profile_sse41_128_64
  #define STATIC static
  #else
  #define PNAME parasail_sw_stats_striped_profile_sse41_128_64
  #define INAME PNAME
  #define STATIC
  #endif
  #endif
  #endif
  
  parasail_result_t* FNAME(
          const char * const restrict s1, const int s1Len,
          const char * const restrict s2, const int s2Len,
          const int open, const int gap, const parasail_matrix_t *matrix)
  {
      parasail_profile_t *profile = parasail_profile_create_stats_sse_128_64(s1, s1Len, matrix);
      parasail_result_t *result = INAME(profile, s2, s2Len, open, gap);
      parasail_profile_free(profile);
      return result;
  }
  
  STATIC parasail_result_t* PNAME(
          const parasail_profile_t * const restrict profile,
          const char * const restrict s2, const int s2Len,
          const int open, const int gap)
  {
      int32_t i = 0;
      int32_t j = 0;
      int32_t k = 0;
      int32_t end_query = 0;
      int32_t end_ref = 0;
      const int s1Len = profile->s1Len;
      const parasail_matrix_t *matrix = profile->matrix;
      const int32_t segWidth = 2; /* number of values in vector unit */
      const int32_t segLen = (s1Len + segWidth - 1) / segWidth;
      __m128i* const restrict vProfile  = (__m128i*)profile->profile64.score;
      __m128i* const restrict vProfileM = (__m128i*)profile->profile64.matches;
      __m128i* const restrict vProfileS = (__m128i*)profile->profile64.similar;
      __m128i* restrict pvHStore        = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHLoad         = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHMStore       = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHMLoad        = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHSStore       = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHSLoad        = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHLStore       = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHLLoad        = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvEStore        = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvELoad         = parasail_memalign___m128i(16, segLen);
      __m128i* const restrict pvEM      = parasail_memalign___m128i(16, segLen);
      __m128i* const restrict pvES      = parasail_memalign___m128i(16, segLen);
      __m128i* const restrict pvEL      = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHMax          = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHMMax          = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHSMax          = parasail_memalign___m128i(16, segLen);
      __m128i* restrict pvHLMax          = parasail_memalign___m128i(16, segLen);
      __m128i vGapO = _mm_set1_epi64x(open);
      __m128i vGapE = _mm_set1_epi64x(gap);
      __m128i vZero = _mm_setzero_si128();
      __m128i vOne = _mm_set1_epi64x(1);
      int64_t score = NEG_INF;
      int64_t matches = NEG_INF;
      int64_t similar = NEG_INF;
      int64_t length = NEG_INF;
      
      __m128i vMaxH = vZero;
      __m128i vMaxHUnit = vZero;
  #ifdef PARASAIL_TABLE
      parasail_result_t *result = parasail_result_new_table3(segLen*segWidth, s2Len);
  #else
  #ifdef PARASAIL_ROWCOL
      parasail_result_t *result = parasail_result_new_rowcol3(segLen*segWidth, s2Len);
      const int32_t offset = (s1Len - 1) % segLen;
      const int32_t position = (segWidth - 1) - (s1Len - 1) / segLen;
  #else
      parasail_result_t *result = parasail_result_new();
  #endif
  #endif
  
      parasail_memset___m128i(pvHMStore, vZero, segLen);
      parasail_memset___m128i(pvHSStore, vZero, segLen);
      parasail_memset___m128i(pvHLStore, vZero, segLen);
      parasail_memset___m128i(pvEM, vZero, segLen);
      parasail_memset___m128i(pvES, vZero, segLen);
      parasail_memset___m128i(pvEL, vZero, segLen);
      parasail_memset___m128i(pvHStore, vZero, segLen);
      parasail_memset___m128i(pvEStore, _mm_set1_epi64x(-open), segLen);
  
      /* outer loop over database sequence */
      for (j=0; j<s2Len; ++j) {
          __m128i vE;
          __m128i vEM;
          __m128i vES;
          __m128i vEL;
          __m128i vF;
          __m128i vFM;
          __m128i vFS;
          __m128i vFL;
          __m128i vH;
          __m128i vHM;
          __m128i vHS;
          __m128i vHL;
          const __m128i* vP = NULL;
          const __m128i* vPM = NULL;
          const __m128i* vPS = NULL;
          __m128i* pv = NULL;
  
          /* Initialize F value to 0.  Any errors to vH values will be corrected
           * in the Lazy_F loop.  */
          vF = vZero;
          vFM = vZero;
          vFS = vZero;
          vFL = vZero;
  
          /* load final segment of pvHStore and shift left by 2 bytes */
          vH = _mm_slli_si128(pvHStore[segLen - 1], 8);
          vHM = _mm_slli_si128(pvHMStore[segLen - 1], 8);
          vHS = _mm_slli_si128(pvHSStore[segLen - 1], 8);
          vHL = _mm_slli_si128(pvHLStore[segLen - 1], 8);
  
          /* Correct part of the vProfile */
          vP = vProfile + matrix->mapper[(unsigned char)s2[j]] * segLen;
          vPM = vProfileM + matrix->mapper[(unsigned char)s2[j]] * segLen;
          vPS = vProfileS + matrix->mapper[(unsigned char)s2[j]] * segLen;
  
          /* Swap the 2 H buffers. */
          pv = pvHLoad;
          pvHLoad = pvHStore;
          pvHStore = pv;
          pv = pvHMLoad;
          pvHMLoad = pvHMStore;
          pvHMStore = pv;
          pv = pvHSLoad;
          pvHSLoad = pvHSStore;
          pvHSStore = pv;
          pv = pvHLLoad;
          pvHLLoad = pvHLStore;
          pvHLStore = pv;
          pv = pvELoad;
          pvELoad = pvEStore;
          pvEStore = pv;
  
          /* inner loop to process the query sequence */
          for (i=0; i<segLen; ++i) {
              __m128i case1not;
              __m128i case2not;
              __m128i case2;
              __m128i case3;
              __m128i cond_zero;
  
              vH = _mm_add_epi64(vH, _mm_load_si128(vP + i));
              vE = _mm_load_si128(pvELoad + i);
  
              /* determine which direction of length and match to
               * propagate, before vH is finished calculating */
              case1not = _mm_or_si128(
                      _mm_cmplt_epi64_rpl(vH,vF),_mm_cmplt_epi64_rpl(vH,vE));
              case2not = _mm_cmplt_epi64_rpl(vF,vE);
              case2 = _mm_andnot_si128(case2not,case1not);
              case3 = _mm_and_si128(case1not,case2not);
  
              /* Get max from vH, vE and vF. */
              vH = _mm_max_epi64_rpl(vH, vE);
              vH = _mm_max_epi64_rpl(vH, vF);
              vH = _mm_max_epi64_rpl(vH, vZero);
              /* Save vH values. */
              _mm_store_si128(pvHStore + i, vH);
              cond_zero = _mm_cmpeq_epi64(vH, vZero);
  
              /* calculate vM */
              vEM = _mm_load_si128(pvEM + i);
              vHM = _mm_blendv_epi8(
                      _mm_add_epi64(vHM, _mm_load_si128(vPM + i)),
                      _mm_or_si128(
                          _mm_and_si128(case2, vFM),
                          _mm_and_si128(case3, vEM)),
                      case1not);
              vHM = _mm_andnot_si128(cond_zero, vHM);
              _mm_store_si128(pvHMStore + i, vHM);
  
              /* calculate vS */
              vES = _mm_load_si128(pvES + i);
              vHS = _mm_blendv_epi8(
                      _mm_add_epi64(vHS, _mm_load_si128(vPS + i)),
                      _mm_or_si128(
                          _mm_and_si128(case2, vFS),
                          _mm_and_si128(case3, vES)),
                      case1not);
              vHS = _mm_andnot_si128(cond_zero, vHS);
              _mm_store_si128(pvHSStore + i, vHS);
  
              /* calculate vL */
              vEL = _mm_load_si128(pvEL + i);
              vHL = _mm_blendv_epi8(
                      _mm_add_epi64(vHL, vOne),
                      _mm_or_si128(
                          _mm_and_si128(case2, _mm_add_epi64(vFL, vOne)),
                          _mm_and_si128(case3, _mm_add_epi64(vEL, vOne))),
                      case1not);
              vHL = _mm_andnot_si128(cond_zero, vHL);
              _mm_store_si128(pvHLStore + i, vHL);
              
  #ifdef PARASAIL_TABLE
              arr_store_si128(result->matches_table, vHM, i, segLen, j, s2Len);
              arr_store_si128(result->similar_table, vHS, i, segLen, j, s2Len);
              arr_store_si128(result->length_table, vHL, i, segLen, j, s2Len);
              arr_store_si128(result->score_table, vH, i, segLen, j, s2Len);
  #endif
              vMaxH = _mm_max_epi64_rpl(vH, vMaxH);
  
              /* Update vE value. */
              vH = _mm_sub_epi64(vH, vGapO);
              vE = _mm_sub_epi64(vE, vGapE);
              vE = _mm_max_epi64_rpl(vE, vH);
              _mm_store_si128(pvEStore + i, vE);
              _mm_store_si128(pvEM + i, vHM);
              _mm_store_si128(pvES + i, vHS);
              _mm_store_si128(pvEL + i, vHL);
  
              /* Update vF value. */
              vF = _mm_sub_epi64(vF, vGapE);
              vF = _mm_max_epi64_rpl(vF, vH);
              vFM = vHM;
              vFS = vHS;
              vFL = vHL;
  
              /* Load the next vH. */
              vH = _mm_load_si128(pvHLoad + i);
              vHM = _mm_load_si128(pvHMLoad + i);
              vHS = _mm_load_si128(pvHSLoad + i);
              vHL = _mm_load_si128(pvHLLoad + i);
          }
  
          /* Lazy_F loop: has been revised to disallow adjecent insertion and
           * then deletion, so don't update E(i, i), learn from SWPS3 */
          for (k=0; k<segWidth; ++k) {
              __m128i vHp = _mm_slli_si128(pvHLoad[segLen - 1], 8);
              vF = _mm_slli_si128(vF, 8);
              vFM = _mm_slli_si128(vFM, 8);
              vFS = _mm_slli_si128(vFS, 8);
              vFL = _mm_slli_si128(vFL, 8);
              for (i=0; i<segLen; ++i) {
                  __m128i case1not;
                  __m128i case2not;
                  __m128i case2;
                  __m128i cond_zero;
  
                  /* need to know where match and length come from so
                   * recompute the cases as in the main loop */
                  vHp = _mm_add_epi64(vHp, _mm_load_si128(vP + i));
                  vE = _mm_load_si128(pvELoad + i);
                  case1not = _mm_or_si128(
                          _mm_cmplt_epi64_rpl(vHp,vF),_mm_cmplt_epi64_rpl(vHp,vE));
                  case2not = _mm_cmplt_epi64_rpl(vF,vE);
                  case2 = _mm_andnot_si128(case2not,case1not);
  
                  vH = _mm_load_si128(pvHStore + i);
                  vH = _mm_max_epi64_rpl(vH,vF);
                  _mm_store_si128(pvHStore + i, vH);
                  cond_zero = _mm_cmpeq_epi64(vH, vZero);
  
                  vHM = _mm_load_si128(pvHMStore + i);
                  vHM = _mm_blendv_epi8(vHM, vFM, case2);
                  vHM = _mm_andnot_si128(cond_zero, vHM);
                  _mm_store_si128(pvHMStore + i, vHM);
                  _mm_store_si128(pvEM + i, vHM);
  
                  vHS = _mm_load_si128(pvHSStore + i);
                  vHS = _mm_blendv_epi8(vHS, vFS, case2);
                  vHS = _mm_andnot_si128(cond_zero, vHS);
                  _mm_store_si128(pvHSStore + i, vHS);
                  _mm_store_si128(pvES + i, vHS);
  
                  vHL = _mm_load_si128(pvHLStore + i);
                  vHL = _mm_blendv_epi8(vHL, _mm_add_epi64(vFL,vOne), case2);
                  vHL = _mm_andnot_si128(cond_zero, vHL);
                  _mm_store_si128(pvHLStore + i, vHL);
                  _mm_store_si128(pvEL + i, vHL);
  
  #ifdef PARASAIL_TABLE
                  arr_store_si128(result->matches_table, vHM, i, segLen, j, s2Len);
                  arr_store_si128(result->similar_table, vHS, i, segLen, j, s2Len);
                  arr_store_si128(result->length_table, vHL, i, segLen, j, s2Len);
                  arr_store_si128(result->score_table, vH, i, segLen, j, s2Len);
  #endif
                  vMaxH = _mm_max_epi64_rpl(vH, vMaxH);
                  vH = _mm_sub_epi64(vH, vGapO);
                  vF = _mm_sub_epi64(vF, vGapE);
                  if (! _mm_movemask_epi8(_mm_cmpgt_epi64_rpl(vF, vH))) goto end;
                  /*vF = _mm_max_epi64_rpl(vF, vH);*/
                  vFM = vHM;
                  vFS = vHS;
                  vFL = vHL;
                  vHp = _mm_load_si128(pvHLoad + i);
              }
          }
  end:
          {
          }
  
          {
              __m128i vCompare = _mm_cmpgt_epi64_rpl(vMaxH, vMaxHUnit);
              if (_mm_movemask_epi8(vCompare)) {
                  score = _mm_hmax_epi64_rpl(vMaxH);
                  vMaxHUnit = _mm_set1_epi64x(score);
                  end_ref = j;
                  (void)memcpy(pvHMax, pvHStore, sizeof(__m128i)*segLen);
                  (void)memcpy(pvHMMax, pvHMStore, sizeof(__m128i)*segLen);
                  (void)memcpy(pvHSMax, pvHSStore, sizeof(__m128i)*segLen);
                  (void)memcpy(pvHLMax, pvHLStore, sizeof(__m128i)*segLen);
              }
          }
  
  #ifdef PARASAIL_ROWCOL
          /* extract last value from the column */
          {
              vH = _mm_load_si128(pvHStore + offset);
              vHM = _mm_load_si128(pvHMStore + offset);
              vHS = _mm_load_si128(pvHSStore + offset);
              vHL = _mm_load_si128(pvHLStore + offset);
              for (k=0; k<position; ++k) {
                  vH = _mm_slli_si128(vH, 8);
                  vHM = _mm_slli_si128(vHM, 8);
                  vHS = _mm_slli_si128(vHS, 8);
                  vHL = _mm_slli_si128(vHL, 8);
              }
              result->score_row[j] = (int64_t) _mm_extract_epi64 (vH, 1);
              result->matches_row[j] = (int64_t) _mm_extract_epi64 (vHM, 1);
              result->similar_row[j] = (int64_t) _mm_extract_epi64 (vHS, 1);
              result->length_row[j] = (int64_t) _mm_extract_epi64 (vHL, 1);
          }
  #endif
      }
  
      /* Trace the alignment ending position on read. */
      {
          int64_t *t = (int64_t*)pvHMax;
          int64_t *m = (int64_t*)pvHMMax;
          int64_t *s = (int64_t*)pvHSMax;
          int64_t *l = (int64_t*)pvHLMax;
          int32_t column_len = segLen * segWidth;
          end_query = s1Len;
          for (i = 0; i<column_len; ++i, ++t, ++m, ++s, ++l) {
              if (*t == score) {
                  int32_t temp = i / segWidth + i % segWidth * segLen;
                  if (temp < end_query) {
                      end_query = temp;
                      matches = *m;
                      similar = *s;
                      length = *l;
                  }
              }
          }
      }
  
  #ifdef PARASAIL_ROWCOL
      for (i=0; i<segLen; ++i) {
          __m128i vH = _mm_load_si128(pvHStore+i);
          __m128i vHM = _mm_load_si128(pvHMStore+i);
          __m128i vHS = _mm_load_si128(pvHSStore+i);
          __m128i vHL = _mm_load_si128(pvHLStore+i);
          arr_store_col(result->score_col, vH, i, segLen);
          arr_store_col(result->matches_col, vHM, i, segLen);
          arr_store_col(result->similar_col, vHS, i, segLen);
          arr_store_col(result->length_col, vHL, i, segLen);
      }
  #endif
  
      
  
      result->score = score;
      result->matches = matches;
      result->similar = similar;
      result->length = length;
      result->end_query = end_query;
      result->end_ref = end_ref;
  
      parasail_free(pvHLMax);
      parasail_free(pvHSMax);
      parasail_free(pvHMMax);
      parasail_free(pvHMax);
      parasail_free(pvEL);
      parasail_free(pvES);
      parasail_free(pvEM);
      parasail_free(pvELoad);
      parasail_free(pvEStore);
      parasail_free(pvHLLoad);
      parasail_free(pvHLStore);
      parasail_free(pvHSLoad);
      parasail_free(pvHSStore);
      parasail_free(pvHMLoad);
      parasail_free(pvHMStore);
      parasail_free(pvHLoad);
      parasail_free(pvHStore);
  
      return result;
  }
  
  #ifdef FASTSTATS
  #ifdef PARASAIL_TABLE
  #else
  #ifdef PARASAIL_ROWCOL
  #else
  #include <assert.h>
  parasail_result_t* INAME(
          const parasail_profile_t * const restrict profile,
          const char * const restrict s2, const int s2Len,
          const int open, const int gap)
  {
      const char *s1 = profile->s1;
      const parasail_matrix_t *matrix = profile->matrix;
  
      /* find the end loc first with the faster implementation */
      parasail_result_t *result = parasail_sw_striped_profile_sse41_128_64(profile, s2, s2Len, open, gap);
      if (!result->saturated) {
  #if 0
          int s1Len_new = 0;
          int s2Len_new = 0;
          char *s1_new = NULL;
          char *s2_new = NULL;
          parasail_profile_t *profile_new = NULL;
          parasail_result_t *result_new = NULL;
          int s1_begin = 0;
          int s2_begin = 0;
          int s1Len_final = 0;
          int s2Len_final = 0;
          parasail_profile_t *profile_final = NULL;
          parasail_result_t *result_final = NULL;
  
          /* using the end loc and the non-stats version of the function,
           * reverse the inputs and find the beg loc */
          s1Len_new = result->end_query+1;
          s2Len_new = result->end_ref+1;
          s1_new = parasail_reverse(s1, s1Len_new);
          s2_new = parasail_reverse(s2, s2Len_new);
          profile_new = parasail_profile_create_sse_128_64(
                  s1_new, s1Len_new, matrix);
          profile_new->stop = result->score;
          result_new = parasail_sw_striped_profile_sse41_128_64(
                  profile_new, s2_new, s2Len_new, open, gap);
  
          /* using both the beg and end loc, call the original stats func */
          s1_begin = s1Len_new - result_new->end_query - 1;
          s2_begin = s2Len_new - result_new->end_ref - 1;
          s1Len_final = s1Len_new - s1_begin;
          s2Len_final = s2Len_new - s2_begin;
          assert(s1_begin >= 0);
          assert(s2_begin >= 0);
          assert(s1Len_new > s1_begin);
          assert(s2Len_new > s2_begin);
          profile_final = parasail_profile_create_stats_sse_128_64(
                  &s1[s1_begin], s1Len_final, matrix);
          result_final = PNAME(
                  profile_final, &s2[s2_begin], s2Len_final, open, gap);
  
          /* clean up all the temporary profiles, sequences, and results */
          free(s1_new);
          free(s2_new);
          parasail_profile_free(profile_new);
          parasail_profile_free(profile_final);
          parasail_result_free(result);
          parasail_result_free(result_new);
  
          /* correct the end locations before returning */
          result_final->end_query = s1Len_new-1;
          result_final->end_ref = s2Len_new-1;
          return result_final;
  #else
          int s1Len_new = 0;
          int s2Len_new = 0;
          parasail_profile_t *profile_final = NULL;
          parasail_result_t *result_final = NULL;
  
          /* using the end loc, call the original stats function */
          s1Len_new = result->end_query+1;
          s2Len_new = result->end_ref+1;
          profile_final = parasail_profile_create_stats_sse_128_64(
                  s1, s1Len_new, matrix);
          result_final = PNAME(
                  profile_final, s2, s2Len_new, open, gap);
  
          /* clean up all the temporary profiles, sequences, and results */
          parasail_profile_free(profile_final);
          parasail_result_free(result);
  
          /* correct the end locations before returning */
          result_final->end_query = s1Len_new-1;
          result_final->end_ref = s2Len_new-1;
          return result_final;
  #endif
      }
      else {
          return result;
      }
  }
  #endif
  #endif
  #endif
  
  
  