From: Heng Li <lh3@live.co.uk>
Date: Sun, 8 May 2011 03:59:17 +0000 (-0400)
Subject: use integer gmax; code cleanup
X-Git-Tag: ksprintf-final~30
X-Git-Url: http://www.kaiwu.me/postgresql/commit/static/gitweb.js?a=commitdiff_plain;h=3f2b354fcc4d4fc00bbf6f90a465983e4e1b5f04;p=klib.git

use integer gmax; code cleanup
---

diff --git a/ksw.c b/ksw.c
index 510f456..aa0e160 100644
--- 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]);