use integer gmax; code cleanup

diff --git a/ksw.c b/ksw.c
index 510f456d3e6b693fd7f1f404926bfdf51da1d959..aa0e160a48cc942449ff310863548bb5a5e92021 100644 (file)
--- a/ksw.c
+++ b/ksw.c
@@ -80,9 +80,9 @@ ksw_query_t *ksw_qinit(int p, int qlen, const uint8_t *query, int m, const int8_
 
 int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) // the first gap costs -(_o+_e)
 {
-       int slen, i, sum, score, m_b, n_b, te = -1;
+       int slen, i, sum, m_b, n_b, te = -1, gmax = 0, _thres = 254 - q->mdiff;
        uint64_t *b;
-       __m128i zero, gapoe, gape, shift, gmax, reduce, thres, *H0, *H1, *E, *Hmax;
+       __m128i zero, gapoe, gape, shift, reduce, thres, *H0, *H1, *E, *Hmax;
 
 #define __max_16(ret, xx) do { \
                (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 8)); \
@@ -94,11 +94,11 @@ int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) /
 
        // initialization
        m_b = n_b = 0; b = 0;
-       gmax = zero = _mm_set1_epi32(0);
+       zero = _mm_set1_epi32(0);
        gapoe = _mm_set1_epi8(a->gapo + a->gape);
        gape = _mm_set1_epi8(a->gape);
        shift = _mm_set1_epi8(q->shift);
-       thres = _mm_set1_epi8(254 - q->mdiff); // when max score exceeds this, shift all scores by "reduce" below
+       thres = _mm_set1_epi8(_thres); // when max score exceeds this, shift all scores by "reduce" below
        reduce = _mm_set1_epi8(127);
        H0 = q->H0; H1 = q->H1; E = q->E; Hmax = q->Hmax;
        q->T = a->T; // avoid potential cache miss
@@ -113,7 +113,7 @@ int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) /
                int j, k, cmp, imax, T = q->T;
                __m128i e, h, f = zero, max = zero, *S = q->qp + target[i] * slen; // s is the 1st score vector
                h = _mm_load_si128(H0 + slen - 1); // h={2,5,8,11,14,17,-1,-1} in the above example
-               h = _mm_slli_si128(h, 1); // h=H(i-1,-1); << instead of >> because x86 is little-endian
+               h = _mm_slli_si128(h, 1); // h=H(i-1,-1); << instead of >> because x64 is little-endian
                for (j = 0; LIKELY(j < slen); ++j) {
                        /* SW cells are computed in the following order:
                         *   H(i,j)   = max{H(i-1,j-1)+S(i,j), E(i,j), F(i,j)}
@@ -153,8 +153,9 @@ int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) /
                        }
                }
 end_loop:
-               __max_16(imax, max); // imax is the maximum number if max
-               if (imax + sum >= T) { // write the b array
+               //int k;for (k=0;k<16;++k)printf("%d ", ((uint8_t*)&max)[k]);printf("\n");
+               __max_16(imax, max); // imax is the maximum number in max
+               if (imax + sum >= T) { // write the b array; this condition adds branching unfornately
                        if (n_b == 0 || (uint32_t)b[n_b-1] + 1 != (uint32_t)i) { // then append
                                if (n_b == m_b) {
                                        m_b = m_b? m_b<<1 : 8;
@@ -162,23 +163,16 @@ end_loop:
                                }
                                b[n_b++] = (uint64_t)(imax + sum)<<32 | i;
                        } else if (b[n_b-1]>>32 < imax + sum) b[n_b-1] = (uint64_t)(imax + sum)<<32 | i; // modify the last
-                       if (_mm_movemask_epi8(_mm_cmpeq_epi8(_mm_subs_epu8(max, gmax), zero)) != 0xffff) { // if max > gmax at any position
-                               int igmax;
-                               __max_16(igmax, gmax); // maximum in gmax
-                               if (imax > igmax) { // if the max has been updated
-                                       te = i; // the end position on the target
-                                       for (j = 0; LIKELY(j < slen); ++j) // keep the H1 vector
-                                               _mm_store_si128(Hmax + j, _mm_load_si128(H1 + j));
-                               }
-                       }
                }
-               gmax = _mm_max_epu8(gmax, max); // update gmax
-               if (_mm_movemask_epi8(_mm_cmpeq_epi8(_mm_subs_epu8(gmax, thres), zero)) != 0xffff) { // reduce
+               if (imax > gmax) {
+                       gmax = imax; te = i; // te is the end position on the target
+                       for (j = 0; LIKELY(j < slen); ++j) // keep the H1 vector
+                               _mm_store_si128(Hmax + j, _mm_load_si128(H1 + j));
+               }
+               if (gmax >= _thres) { // reduce
                        // When overflow is going to happen, subtract 127 from all scores. This heuristic
                        // may miss the best alignment, but in practice this should happen very rarely.
-                       sum += 127;
-                       gmax = _mm_subs_epu8(gmax, reduce);
-                       //int k;for (k=0;k<16;++k)printf("%d ", ((uint8_t*)&thres)[k]);printf("\n");
+                       sum += 127; gmax -= 127;
                        for (j = 0; LIKELY(j < slen); ++j) {
                                h = _mm_subs_epu8(_mm_load_si128(H1 + j), reduce);
                                _mm_store_si128(H1 + j, h);
@@ -188,16 +182,15 @@ end_loop:
                }
                S = H1; H1 = H0; H0 = S; // swap H0 and H1
        }
-       __max_16(score, gmax);
-       a->score = score + sum; a->te = te;
+       a->score = gmax + sum; a->te = te;
        { // get a->qe, the end of query match
                int max = -1;
                uint8_t *t = (uint8_t*)Hmax;
-               for (i = 0, score = 0, a->qe = -1; i < q->qlen; ++i, ++t)
+               for (i = 0, a->qe = -1; i < q->qlen; ++i, ++t)
                        if ((int)*t > max) max = *t, a->qe = i / 16 + i % 16 * slen;
        }
        free(b);
-       return a->score < a->T? -1 : a->score;
+       return a->score;
 }
 
 /*******************************************
@@ -284,10 +277,10 @@ int main(int argc, char *argv[])
                        int s;
                        for (i = 0; i < kst->seq.l; ++i) kst->seq.s[i] = seq_nt4_table[(int)kst->seq.s[i]];
                        s = ksw_sse2_16(q[0], kst->seq.l, (uint8_t*)kst->seq.s, &a);
-                       if (s >= 0) printf("%s\t%s\t+\t%d\t%d\t%d\n", ksq->name.s, kst->name.s, s, a.te+1, a.qe+1);
+                       printf("%s\t%s\t+\t%d\t%d\t%d\n", ksq->name.s, kst->name.s, s, a.te+1, a.qe+1);
                        if (q[1]) {
                                s = ksw_sse2_16(q[1], kst->seq.l, (uint8_t*)kst->seq.s, &a);
-                               if (s >= 0) printf("%s\t%s\t-\t%d\t%d\t%d\n", ksq->name.s, kst->name.s, s, a.te+1, a.qe+1);
+                               printf("%s\t%s\t-\t%d\t%d\t%d\n", ksq->name.s, kst->name.s, s, a.te+1, a.qe+1);
                        }
                }
                free(q[0]); free(q[1]);