From ba2f19ba377d9337a742bc7bab9b89e509d2a1e4 Mon Sep 17 00:00:00 2001 From: Saurabh Date: Thu, 18 Nov 2021 21:14:38 -0800 Subject: [PATCH] mm2-fast-v2.22 init --- Makefile | 42 +- align.c | 66 ++ build_rmi.sh | 16 + index.c | 360 +++++++ ksw2_extd2_avx.c | 2319 ++++++++++++++++++++++++++++++++++++++++++++++ ksw2_extd2_avx.h | 42 + lchain.c | 218 ++++- lisa_hash.h | 307 ++++++ main.c | 100 +- map.c | 20 + minimap.h | 20 + seed.c | 45 + 12 files changed, 3535 insertions(+), 20 deletions(-) create mode 100755 build_rmi.sh create mode 100644 ksw2_extd2_avx.c create mode 100644 ksw2_extd2_avx.h create mode 100644 lisa_hash.h diff --git a/Makefile b/Makefile index 4118616..1d0fc3f 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,42 @@ CFLAGS= -g -Wall -O2 -Wc++-compat #-Wextra -CPPFLAGS= -DHAVE_KALLOC -INCLUDES= +CPPFLAGS= -DHAVE_KALLOC #-march=native #-DALIGN_AVX -DPARALLEL_CHAINING #-DMANUAL_PROFILING +COMP_FLAG = -march=native + +ifeq ($(avx2_compile), 1) + COMP_FLAG = -mavx2 +endif + +#CPPFLAGS= -DHAVE_KALLOC -mavx2 -DALIGN_AVX -DAPPLY_AVX2 -DPARALLEL_CHAINING #-DLISA_HASH -DUINT64 -DVECTORIZE #-DMANUAL_PROFILING +#CPPFLAGS= -DHAVE_KALLOC -mavx2 -DPARALLEL_CHAINING #-DMANUAL_PROFILING + +OPT_FLAGS= -DPARALLEL_CHAINING -DALIGN_AVX -DAPPLY_AVX2 +OPT_FLAGS+=$(COMP_FLAG) +ifeq ($(lhash_index), 1) + CPPFLAGS+= -DLISA_INDEX +endif +ifeq ($(lhash), 1) + OPT_FLAGS+= -DLISA_HASH -DUINT64 -DVECTORIZE +endif +ifeq ($(manual_profile), 1) + CPPFLAGS+= -DMANUAL_PROFILING +endif + +#ifeq ($(use_avx2), 1) +# OPT_FLAGS+= -DAPPLY_AVX2 +#endif + +ifeq ($(disable_output), 1) + CPPFLAGS+= -DDISABLE_OUTPUT +endif + +ifeq ($(no_opt),) + CPPFLAGS+= $(OPT_FLAGS) +endif + + + +#INCLUDES= +#INCLUDES= -I./ext/TAL_offline/src/LISA-hash #-I./ext/TAL/src/dynamic-programming OBJS= kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o \ lchain.o align.o hit.o seed.o map.o format.o pe.o esterr.o splitidx.o \ ksw2_ll_sse.o @@ -10,7 +46,7 @@ LIBS= -lm -lz -lpthread ifeq ($(arm_neon),) # if arm_neon is not defined ifeq ($(sse2only),) # if sse2only is not defined - OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o + OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o ksw2_extd2_avx.o else # if sse2only is defined OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o endif diff --git a/align.c b/align.c index ddbb0bd..d0e3a10 100644 --- a/align.c +++ b/align.c @@ -5,6 +5,13 @@ #include "minimap.h" #include "mmpriv.h" #include "ksw2.h" +#include "ksw2_extd2_avx.h" +#include +extern uint64_t avg; +extern uint64_t alignment_time; +extern void *km1; +extern uint64_t km_size;// = 500000000; // 500 MB +extern int km_top; static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc_ambi) { @@ -313,6 +320,7 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // } } +#if 0 static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez) { if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) { @@ -340,7 +348,65 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint fprintf(stderr, "\n"); } } +#endif +#if 1 +static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez) +{ +#ifdef MANUAL_PROFILING + uint64_t align_start = __rdtsc(); +#endif + if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) { + int i; + fprintf(stderr, "===> q=(%d,%d), e=(%d,%d), bw=%d, flag=%d, zdrop=%d <===\n", opt->q, opt->q2, opt->e, opt->e2, w, flag, opt->zdrop); + for (i = 0; i < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr); + fputc('\n', stderr); + for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr); + fputc('\n', stderr); + } + if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) { + ksw_reset_extz(ez); + ez->zdropped = 1; + } else if (opt->flag & MM_F_SPLICE) + ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, opt->junc_bonus, flag, junc, ez); + else if (opt->q == opt->q2 && opt->e == opt->e2) + ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez); + else{ +#if defined (ALIGN_AVX) && (defined(__AVX512BW__) || (defined(__AVX2__) && defined(APPLY_AVX2))) +#ifdef __AVX512BW__ + + ksw_extd2_avx512(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez); +#elif __AVX2__ + avg = 0; +// uint64_t *ptr_km = (uint64_t *) km1; +// for(uint64_t itr = 0; itr < km_size/512; itr++){ +// avg+=ptr_km[itr]; +// } + +//#ifdef MANUAL_PROFILING +// uint64_t align_start = __rdtsc(); +//#endif + ksw_extd2_avx2(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez); +//#ifdef MANUAL_PROFILING +// alignment_time += (__rdtsc() - align_start); +//#endif +#endif +#else + ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez); +#endif + } + if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) { + int i; + fprintf(stderr, "score=%d, cigar=", ez->score); + for (i = 0; i < ez->n_cigar; ++i) + fprintf(stderr, "%d%c", ez->cigar[i]>>4, "MIDN"[ez->cigar[i]&0xf]); + fprintf(stderr, "\n"); + } +#ifdef MANUAL_PROFILING + alignment_time += (__rdtsc() - align_start); +#endif +} +#endif static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x) { int64_t i, off0 = mi->seq[rid].offset, off = off0 + x; diff --git a/build_rmi.sh b/build_rmi.sh new file mode 100755 index 0000000..396dbd4 --- /dev/null +++ b/build_rmi.sh @@ -0,0 +1,16 @@ +ref_data=$1 +preset=$2 + +make clean && make lhash_index=1 +touch temp_read.fastq +./minimap2 -ax $2 $1 temp_read.fastq >/dev/null + +kv_file=$1"_"$2"_minimizers_key_value_sorted" + +full_path=`readlink -f $kv_file` + +cd ./ext/TAL_offline +make lisa_hash +./build-lisa-hash-index $full_path + +rm ../../temp_read.fastq diff --git a/index.c b/index.c index 0d8a2ae..6e06e77 100644 --- a/index.c +++ b/index.c @@ -14,6 +14,21 @@ #include "mmpriv.h" #include "kvec.h" #include "khash.h" +#include +#include +#include +#include +#include +using namespace std; + + + +#ifdef LISA_HASH +#include "lisa_hash.h" +extern lisa_hash *lh; +#endif + + #define idx_hash(a) ((a)>>1) #define idx_eq(a, b) ((a)>>1 == (b)>>1) @@ -52,9 +67,43 @@ mm_idx_t *mm_idx_init(int w, int k, int b, int flag) if (!(mm_dbg_flag & 1)) mi->km = km_init(); return mi; } +void mm_idx_destroy_mm_hash(mm_idx_t *mi) +{ + //fprintf(stderr, "mm_destroy_hash\n"); + uint32_t i; + if (mi == 0) return; + if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h); + if (mi->B) { + for (i = 0; i < 1U<b; ++i) { + free(mi->B[i].p); + free(mi->B[i].a.a); + kh_destroy(idx, (idxhash_t*)mi->B[i].h); + } + } +} +void mm_idx_destroy_seq(mm_idx_t *mi) +{ + //fprintf(stderr, "mm_destroy_seq\n"); + + uint32_t i; + if (mi == 0) return; + if (mi->I) { + for (i = 0; i < mi->n_seq; ++i) + free(mi->I[i].a); + free(mi->I); + } + if (!mi->km) { + for (i = 0; i < mi->n_seq; ++i) + free(mi->seq[i].name); + free(mi->seq); + } else km_destroy(mi->km); + free(mi->B); free(mi->S); free(mi); +} + void mm_idx_destroy(mm_idx_t *mi) { + uint32_t i; if (mi == 0) return; if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h); @@ -96,6 +145,317 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n) return &b->p[kh_val(h, k)>>32]; } } +//Output minimap2's hash table entries +class hash_entry { + public: + uint64_t key; + uint64_t n; + uint64_t *p; + hash_entry(uint64_t k, uint64_t n_, uint64_t *p_){ + key = k; + n = n_; + p = p_; + } + +}; +bool key_sort( hash_entry i1, hash_entry i2) +{ + return (i1.key < i2.key); +} + +#if 0 +void mm_idx_load_key_value_lisa(const char* f_name, const mm_idx_t *mi) +{ + uint64_t tic = __rdtsc(); + std::vector v_hash; + + //ofstream f(f_name); + fprintf(stderr, "Building sorted key-val map\n"); + + uint32_t i,j; + uint64_t num_values = 0; + for (i = 0; i < 1U<b; ++i) { + + + //fprintf(stderr, "BucketID %lu \n", i); + idxhash_t *h = (idxhash_t*)mi->B[i].h; + khint_t k; + if (h == 0) continue; + for (k = 0; k < kh_end(h); ++k){ + if (kh_exist(h, k)) { + uint64_t key = kh_key(h, k), bucket_id = i; + key = key>>1; + + key = key<b | bucket_id; + + if(kh_key(h, k)&1) + { + //print key value + //fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0); + v_hash.push_back(hash_entry(key, kh_val(h, k), NULL)); + } + else + { // print key + uint32_t n = (uint32_t)kh_val(h, k); + //fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n); + // for 0 to lsb 32 val + // print b->p[msb 32 of val] + + v_hash.push_back(hash_entry(key, n, &mi->B[i].p[(kh_val(h, k)>>32) + 0])); + } + } + } + } + sort(v_hash.begin(), v_hash.end(), key_sort); + fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic); + fprintf(stderr, "Storing hash to %s \n", f_name); + tic = __rdtsc(); + + int64_t itr_p = 0; + for( int i = 0; i < v_hash.size(); i++){ + + if(v_hash[i].p == NULL){ + //f<p[itr_p++] = v_hash[i].n; + continue; + } + + //f<p[itr_p++] = v_hash[i].p[j]; + num_values++; + } + //f<> m; + std::vector v_hash; + + //ofstream f(f_name); + fprintf(stderr, "Building sorted key-val map\n"); + + uint32_t i,j; + uint64_t num_values = 0; + for (i = 0; i < 1U<b; ++i) { + + + //fprintf(stderr, "BucketID %lu \n", i); + idxhash_t *h = (idxhash_t*)mi->B[i].h; + khint_t k; + if (h == 0) continue; + for (k = 0; k < kh_end(h); ++k){ + if (kh_exist(h, k)) { + uint64_t key = kh_key(h, k), bucket_id = i; + key = key>>1; + + key = key<b | bucket_id; + + if(kh_key(h, k)&1) + { + //print key value + //fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0); + //m[key].push_back(kh_val(h, k)); + v_hash.push_back(hash_entry(key, kh_val(h, k), NULL)); + } + else + { // print key + uint32_t n = (uint32_t)kh_val(h, k); + //fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n); + // for 0 to lsb 32 val + // print b->p[msb 32 of val] + + v_hash.push_back(hash_entry(key, n, &mi->B[i].p[(kh_val(h, k)>>32) + 0])); + } + } + } + } + sort(v_hash.begin(), v_hash.end(), key_sort); + fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic); + fprintf(stderr, "Storing hash to %s \n", f_name); + tic = __rdtsc(); + + vector key_list; + vector val_list; + vector p_list; +/* + key_list.push_back(m.size()); + for(auto k : m){ + key_list.push_back(k.first); + f<> m; + + ofstream f(f_name); + fprintf(stderr, "Building sorted key-val map\n"); + + uint32_t i,j; + uint64_t num_values = 0; + for (i = 0; i < 1U<b; ++i) { + + + //fprintf(stderr, "BucketID %lu \n", i); + idxhash_t *h = (idxhash_t*)mi->B[i].h; + khint_t k; + if (h == 0) continue; + for (k = 0; k < kh_end(h); ++k){ + if (kh_exist(h, k)) { + uint64_t key = kh_key(h, k), bucket_id = i; + key = key>>1; + + key = key<b | bucket_id; + + if(kh_key(h, k)&1) + { + //print key value + //fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0); + m[key].push_back(kh_val(h, k)); + } + else + { // print key + uint32_t n = (uint32_t)kh_val(h, k); + //fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n); + // for 0 to lsb 32 val + // print b->p[msb 32 of val] + for(j = 0; j < n; j++) + { + //fprintf(stderr, "%llu ", mi->B[i].p[(kh_val(h, k)>>32) + j]); + m[key].push_back(mi->B[i].p[(kh_val(h, k)>>32) + j]); + } + } + } + } + } + fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic); + fprintf(stderr, "Storing hash to %s \n", f_name); + tic = __rdtsc(); + vector key_list; + key_list.push_back(m.size()); + for(auto k : m){ + key_list.push_back(k.first); + f<; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ +#include "ksw2_extd2_avx.h" + +#ifdef __AVX512BW__ + #define __dp_code_block1_pcl_512_mm2_fast \ + /*__mmask64 mskc_ = mskc; */ \ + /*if (t == st_) */ \ + __mmask64 mskc_ = (t == st_) ? mskc_ar[(st0 - t*64)/16]:mskc; \ + z = _mm512_load_si512(&s[t]); \ + xt1 = _mm512_load_si512(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \ + /* tmp = _mm_srli_si128(xt1, 15); */ /* tmp <- x[r-1][t+15] */ \ + tmp = _mm512_set1_epi8(((int8_t*)x)[t*64 + 63]); \ + /* xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_);*/ /* xt1 <- x[r-1][t-1..t+14] */ \ + shf512a = _mm512_shuffle_epi8(xt1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + xt1 = _mm512_mask_blend_epi8(mskc_, slli512, x1_); \ + x1_ = tmp; \ + vt1 = _mm512_load_si512(&v[t]); /* vt1 <- v[r-1][t..t+15] */ \ + /* tmp = _mm_srli_si128(vt1, 15); */ /* tmp <- v[r-1][t+15] */ \ + tmp = _mm512_set1_epi8(((int8_t*)v)[t*64 + 63]); \ + /* vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); *//* vt1 <- v[r-1][t-1..t+14] */ \ + shf512a = _mm512_shuffle_epi8(vt1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + vt1 = _mm512_mask_blend_epi8(mskc_, slli512, v1_); \ + v1_ = tmp; \ + a = _mm512_add_epi8(xt1, vt1); /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \ + ut = _mm512_load_si512(&u[t]); /* ut <- u[t..t+15] */ \ + b = _mm512_add_epi8(_mm512_load_si512(&y[t]), ut); /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */ \ + x2t1= _mm512_load_si512(&x2[t]); \ + /* tmp = _mm_srli_si128(x2t1, 15);*/ \ + tmp = _mm512_set1_epi8(((int8_t*)x2)[t*64 + 63]); \ + /* x2t1= _mm_or_si128(_mm_slli_si128(x2t1, 1), x21_); */ \ + shf512a = _mm512_shuffle_epi8(x2t1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + x2t1 = _mm512_mask_blend_epi8(mskc_, slli512, x21_); \ + x21_= tmp; \ + a2= _mm512_add_epi8(x2t1, vt1); \ + b2= _mm512_add_epi8(_mm512_load_si512(&y2[t]), ut); + + +void ksw_extd2_avx512(void *km, int qlen, const uint8_t *query, int tlen, + const uint8_t *target, int8_t m, const int8_t *mat, + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, + int zdrop, int end_bonus, int flag, ksw_extz_t *ez) +{ + + // printf("AVX512 call\n"); + + __m512i ind512_srli = _mm512_set1_epi8(15); + __mmask8 msk_srli = 0x0003; + + + int8_t index[64] __attribute((aligned(64))); + + for (int i=0; i<64; i++) + index[i] = i%16 - 1; + + index[0] = 15; + index[16] = 31; + index[32] = 47; + index[48] = 63; + + __m512i shf512a, shf512b, slli512; + __m512i ind512_slli = _mm512_load_si512((__m512i*) index); + __mmask8 mska = 0x90; + __mmask64 mskb = 0x0001000100010000; + __mmask64 mskc = 0x1; + __mmask64 mskc_ar[4] = {0x1, 0x10000, 0x100000000, 0x1000000000000}; + + + + + __mmask64 msk_ar[5] = {0xFFFF, 0xFFFFFFFF, 0xFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF}; + __mmask64 msk_ar2[5] = {0xFFFF, 0xFFFF, 0xFFFFFFFF, 0xFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF}; + + int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc, long_thres, long_diff; + int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX); + int32_t *H = 0, H0 = 0, last_H0_t = 0; + uint8_t *qr, *sf, *mem, *mem2 = 0; + + __m512i q_, q2_, qe_, qe2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_; + __m512i *u, *v, *x, *y, *x2, *y2, *s, *p = 0; + __m512i one_, two_, three_, four_, s1_, s2_, s3_, s4_; + + ksw_reset_extz(ez); + if (m <= 1 || qlen <= 0 || tlen <= 0) return; + + if (q2 + e2 < q + e) t = q, q = q2, q2 = t, t = e, e = e2, e2 = t; // make sure q+e no larger than q2+e2 + s1_ = _mm512_set1_epi8(0x08); + s2_ = _mm512_set1_epi8(0x10); + s3_ = _mm512_set1_epi8(0x20); + s4_ = _mm512_set1_epi8(0x40); + + one_ = _mm512_set1_epi8(1); + two_ = _mm512_set1_epi8(2); + three_ = _mm512_set1_epi8(3); + four_ = _mm512_set1_epi8(4); + __m512i four32_ = _mm512_set1_epi32(4); + __m512i bt32_ = _mm512_setr_epi32(0,0,0,0,4,4,4,4,8,8,8,8,12,12,12,12); + + zero_ = _mm512_set1_epi8(0); + q_ = _mm512_set1_epi8(q); + q2_ = _mm512_set1_epi8(q2); + qe_ = _mm512_set1_epi8(q + e); + qe2_ = _mm512_set1_epi8(q2 + e2); + sc_mch_ = _mm512_set1_epi8(mat[0]); + sc_mis_ = _mm512_set1_epi8(mat[1]); + sc_N_ = mat[m*m-1] == 0? _mm512_set1_epi8(-e2) : _mm512_set1_epi8(mat[m*m-1]); + m1_ = _mm512_set1_epi8(m - 1); // wildcard + + if (w < 0) w = tlen > qlen? tlen : qlen; + wl = wr = w; + tlen_ = (tlen + 63) / 64; + n_col_ = qlen < tlen? qlen : tlen; + n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 63) / 64 + 1; + qlen_ = (qlen + 63) / 64; + for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) { + max_sc = max_sc > mat[t]? max_sc : mat[t]; + min_sc = min_sc < mat[t]? min_sc : mat[t]; + } + if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches + + long_thres = e != e2? (q2 - q) / (e - e2) - 1 : 0; + if (q2 + e2 + long_thres * e2 > q + e + long_thres * e) + ++long_thres; + long_diff = long_thres * (e - e2) - (q2 - q) - e2; + + mem = (uint8_t*)kcalloc(km, tlen_ * 8 + qlen_ + 1 + 63, 64); //TODO: saurabh - 63 + u = (__m512i*)(((size_t)mem + 63) >> 6 << 6); // 16-byte aligned //TODO: saurabh mm_malloc + v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_; + s = y2 + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 64; + memset(u, -q - e, tlen_ * 64); + memset(v, -q - e, tlen_ * 64); + memset(x, -q - e, tlen_ * 64); + memset(y, -q - e, tlen_ * 64); + memset(x2, -q2 - e2, tlen_ * 64); + memset(y2, -q2 - e2, tlen_ * 64); + if (!approx_max) { + H = (int32_t*)kmalloc(km, tlen_ * 64 * 4); + for (t = 0; t < tlen_ * 64; ++t) H[t] = KSW_NEG_INF; + } + if (with_cigar) { + // printf("qlen: %dm tlen: %dm ncol: %d\n", qlen, tlen, n_col_); + mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 64); + p = (__m512i*)(((size_t)mem2 + 63) >> 6 << 6);//TODO: saurabh + off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2); + off_end = off + qlen + tlen - 1; + } + + + #define ORIG 0 // for B1 block optimization + for (t = 0; t < qlen; ++t) { + #if ORIG + qr[t] = query[qlen - 1 - t]; + #else + if (query[qlen - 1 - t] == 4) + qr[t] = 8; + else + qr[t] = query[qlen - 1 - t]; + #endif + //assert(qr[t] < 5); + } + memcpy(sf, target, tlen); + + + #if TIMER_O + uint64_t tic = _rdtsc(); + #endif + + #if 1 // changes due to B1 modification + int8_t cN = mat[m*m-1] == 0? -e2 : mat[m*m-1]; + int8_t pmat[64] __attribute((aligned(64))) = {0}; + pmat[0] = mat[0]; + pmat[1] = pmat[2] = pmat[3] = mat[1]; + pmat[4] = pmat[5] = pmat[6] = pmat[7] = pmat[8] = pmat[9] = pmat[10] = pmat[11] = pmat[12] = cN; + __m512i vpmat = _mm512_load_si512(pmat); + vpmat = _mm512_shuffle_i32x4(vpmat, vpmat, 0x0000); + #endif + + static int itr = 0; + itr++; + + for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) { + int st = 0, en = tlen - 1, st0, en0, st_, en_; + int8_t x1, x21, v1; + uint8_t *qrr = qr + (qlen - 1 - r); + int8_t *u8 = (int8_t*)u, *v8 = (int8_t*)v, *x8 = (int8_t*)x, *x28 = (int8_t*)x2; + __m512i x1_, x21_, v1_; + // find the boundaries + if (st < r - qlen + 1) st = r - qlen + 1; + if (en > r) en = r; + if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil + if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor + if (st > en) { + ez->zdropped = 1; + break; + } + st0 = st, en0 = en; + { + // _mm_prefetch(sf + st0, _MM_HINT_T0); + // _mm_prefetch(qrr + st0, _MM_HINT_T0); + // _mm_prefetch(s + st0, _MM_HINT_T0); + // _mm_prefetch(sf + st0 + 64, _MM_HINT_T0); + // _mm_prefetch(qrr + st0 + 64, _MM_HINT_T0); + // _mm_prefetch(s + st0 + 64, _MM_HINT_T0); + } + + int st_new = st / 16 * 16, en_new = (en + 16) / 16 * 16 - 1;//TODO: saurabh + // int st_new = st / 64 * 64, en_new = (en + 64) / 64 * 64 - 1;//TODO: saurabh + int stb = st, enb = en; + st = st / 64 * 64, en = (en + 64) / 64 * 64 - 1;//TODO: saurabh + int stn = stb / 16 * 16, enn = (enb + 16) / 16 * 16 - 1;//TODO: saurabh + // set boundary conditions + //if(st_new != st) + { + //printf("st0 %d en0 %d st16 %d en16 %d st64 %d en64 %d\n", st0, en0, st_new, en_new, st, en); + //processed_vect64 += ceil(((double)(en - st) + 1)/64); possible_vect64 += ceil(((double)(en_new - st_new) + 1)/64); + //exit(0); + } + if (st_new > 0) { + if (st_new - 1 >= last_st && st_new - 1 <= last_en) { + x1 = x8[st_new - 1], x21 = x28[st_new - 1], v1 = v8[st_new - 1]; // (r-1,s-1) calculated in the last round + } else { + x1 = -q - e, x21 = -q2 - e2; + v1 = -q - e; + } + } else { + x1 = -q - e, x21 = -q2 - e2; + v1 = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; + } + if (en_new >= r) { + ((int8_t*)y)[r] = -q - e, ((int8_t*)y2)[r] = -q2 - e2; + u8[r] = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; + } + #if TIMER_I + uint64_t tic = _rdtsc(); + #endif + + + + // loop fission: set scores first //TODO: saurabh: introducing mask64 + if (!(flag & KSW_EZ_GENERIC_SC)) { + #if ORIG + for (t = st0; t <= en0; t += 64) { + __m512i sq, st, tmp_512, mask_512; + __mmask64 tmp, mask; + sq = _mm512_loadu_si512((__m512i*)&sf[t]); + st = _mm512_loadu_si512((__m512i*)&qrr[t]); + mask = (_mm512_cmpeq_epi8_mask(sq, m1_) | _mm512_cmpeq_epi8_mask(st, m1_)); + tmp = _mm512_cmpeq_epi8_mask(sq, st); + + #if 1 + tmp_512 = _mm512_mask_blend_epi8(tmp, sc_mis_, sc_mch_); + tmp_512 = _mm512_mask_blend_epi8(mask, tmp_512, sc_N_); + #else + tmp = _mm_or_si128(_mm_andnot_si128(tmp, sc_mis_), _mm_and_si128(tmp, sc_mch_));//TODO: saurabh - left unchanged + tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp), _mm_and_si128(mask, sc_N_)); + #endif + if (t + 64 > en0) + { + __mmask64 msk; + int ind = (en0 - t + 16) >> 4; + // assert(ind >= 0 && ind < 5); + msk = msk_ar2[ind]; + _mm512_mask_storeu_epi8((__m512i*)((int8_t*)s + t), msk, tmp_512); + } + else + _mm512_storeu_si512((__m512i*)((int8_t*)s + t), tmp_512); + } + #else + for (t = 0; t <= (en0 - st0) - 64; t += 64) { + __m512i sq, st, tmp_512, mask_512; + __mmask64 tmp, mask; + // const int pdf = 4; + //_mm_prefetch(sf + t + 64*pdf, _MM_HINT_T0); + //_mm_prefetch(qrr + t + 64*pdf, _MM_HINT_T0); + //_mm_prefetch(s + t + 64*pdf, _MM_HINT_T0); + + sq = _mm512_loadu_si512((__m512i*)&sf[st0 + t]); + st = _mm512_loadu_si512((__m512i*)&qrr[st0 + t]); + tmp_512 = _mm512_xor_si512(sq, st); + tmp_512 = _mm512_shuffle_epi8(vpmat, tmp_512); + _mm512_storeu_si512((__m512i*)((int8_t*)s + st0 + t), tmp_512); + } + { + __m512i sq, st, tmp_512, mask_512; + __mmask64 tmp, mask; + sq = _mm512_loadu_si512((__m512i*)&sf[st0 + t]); + st = _mm512_loadu_si512((__m512i*)&qrr[st0 + t]); + tmp_512 = _mm512_xor_si512(sq, st); + tmp_512 = _mm512_shuffle_epi8(vpmat, tmp_512); + + if (t + 64 > en0 - st0) + { + __mmask64 msk; + int ind = ((en0 - st0) - t + 16) >> 4; + msk = msk_ar2[ind]; + _mm512_mask_storeu_epi8((__m512i*)((int8_t*)s + st0 + t), msk, tmp_512); + } + else + _mm512_storeu_si512((__m512i*)((int8_t*)s + st0 + t), tmp_512); + } + #endif + } + else + { + // printf("ELSE\n"); + for (t = st0; t <= en0; ++t) { + #if ORIG + ((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]]; + #else + int8_t c = qrr[t]; + if (qrr[t] == 8) { + c = 4; + } + ((uint8_t*)s)[t] = mat[sf[t] * m + c]; + #endif + } + } + + + #if TIMER_I + uint64_t toc = _rdtsc(); + tprof[B1][0] += toc - tic; + #endif + //continue; + // core loop //TODO: saurabh - need more understanding of logic + // fprintf(stderr, "- r: %d, x1: %d, x21: %d, v1: %d, en_new: %d, e: %d, q: %d\n", + //r, x1, x21, v1, en_new, e, q); + x1_ = _mm512_set1_epi8((uint8_t)x1); + x21_ = _mm512_set1_epi8((uint8_t)x21); + v1_ = _mm512_set1_epi8((uint8_t)v1); + + //st_ = st / 16, en_ = en / 16; + st_ = st / 64, en_ = en / 64; + //assert(en_ - st_ + 1 <= n_col_); + if (!with_cigar) { // score only + #if TIMER_I + uint64_t tic = _rdtsc(); + #endif + for (t = st_; t <= en_; ++t) { + __m512i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + + //__dp_code_block1_pcl; + __dp_code_block1_pcl_512_mm2_fast; + z = _mm512_max_epi8(z, a); + z = _mm512_max_epi8(z, b); + z = _mm512_max_epi8(z, a2); + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; // save u[] and v[]; update a, b, a2 and b2 + #if 1 + if (t == en_) { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + // fprintf(stderr, "en0: %d, t: %d, ind: %d, msk: %d\n", en0, t, ind, msk); + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + + } + else { + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + #endif + + #if 1 + if (t == en_) { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + //assert(ind >= 0); + msk = msk_ar2[ind]; + // fprintf(stderr, "en0: %d, t: %d, ind: %d, msk: %d\n", en0, t, ind, msk); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_max_epi8(a, zero_), qe_)); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_max_epi8(b, zero_), qe_)); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_max_epi8(a2, zero_), qe2_)); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_max_epi8(b2, zero_), qe2_)); + } + else + #endif + { + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_max_epi8(a, zero_), qe_)); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_max_epi8(b, zero_), qe_)); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_max_epi8(a2, zero_), qe2_)); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_max_epi8(b2, zero_), qe2_)); + } + // for (int l=0; l<64; l++) + // fprintf(stderr, "%d ", ((int8_t*)x)[l]); + + } + #if TIMER_I + uint64_t toc = _rdtsc(); + tprof[B2][0] += toc - tic; + #endif + + } else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment + __m512i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + + //_mm_prefetch(pr, _MM_HINT_T0); + for (t = st_; t < en_; ++t) { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + //__dp_code_block1_pcl; + __dp_code_block1_pcl_512_mm2_fast; + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a, z),zero_, one_); // d = a > z? 1 : 0 + z = _mm512_max_epi8(z, a); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b, z), d, two_); // d = b > z? 2 : d + z = _mm512_max_epi8(z, b); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a2, z), d, three_); // d = a2 > z? 3 : d + z = _mm512_max_epi8(z, a2); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b2, z), d, four_); // d = b2 > z? 4 : d + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + tmp_mask = _mm512_cmpgt_epi8_mask(a, zero_); + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s1_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b, zero_); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s2_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(a2, zero_); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s3_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b2, zero_); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s4_)); // d = b > 0? 1<<6 : 0 + _mm512_store_si512(&pr[t], d); + //_mm_prefetch(&pr[t+1], _MM_HINT_T0); + } + #if TIMER_I + uint64_t tic = _rdtsc(); + #endif + { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + //__dp_code_block1_pcl; + __dp_code_block1_pcl_512_mm2_fast; + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a, z),zero_, one_); // d = a > z? 1 : 0 + z = _mm512_max_epi8(z, a); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b, z), d, two_); // d = b > z? 2 : d + z = _mm512_max_epi8(z, b); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a2, z), d, three_); // d = a2 > z? 3 : d + z = _mm512_max_epi8(z, a2); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b2, z), d, four_); // d = b2 > z? 3 : d + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + #if 1 + { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + // assert(ind >= 0 && ind < 5); + msk = msk_ar2[ind]; + + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + #endif + + #if 1 + { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + off_end[r] -= (4-ind)*16; + + tmp_mask = _mm512_cmpgt_epi8_mask(a, zero_); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s1_)); // d = a > 0? 1<<3 : 0 + + tmp_mask = _mm512_cmpgt_epi8_mask(b, zero_); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s2_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(a2, zero_); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s3_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b2, zero_); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s4_)); // d = b > 0? 1<<6 : 0 + //_mm512_store_si512(&pr[t], d); + _mm512_mask_storeu_epi8(&pr[t], msk, d); + } + #endif + } + #if TIMER_I + uint64_t toc = _rdtsc(); + tprof[B3][0]+= toc - tic; + #endif +/* + for (t = st_; t <= en_; ++t) { + + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a, z), _mm512_set1_epi8(0), _mm512_set1_epi8(1)); // d = a > z? 1 : 0 + z = _mm512_max_epi8(z, a); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b, z), d, _mm512_set1_epi8(2)); // d = b > z? 2 : d + z = _mm512_max_epi8(z, b); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a2, z), d, _mm512_set1_epi8(3)); // d = a2 > z? 3 : d + z = _mm512_max_epi8(z, a2); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b2, z), d, _mm512_set1_epi8(4)); // d = b2 > z? 3 : d + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + #if 1 + if (t == en_) { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + //assert(ind >= 0 && ind < 5); + msk = msk_ar2[ind]; + + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + + } + else { + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + #endif + + #if 1 + if (t == en_) { + __mmask64 msk; + int ind = (en0 - t*64 + 16) >> 4; + msk = msk_ar2[ind]; + off_end[r] -= (4-ind)*16; + + tmp_mask = _mm512_cmpgt_epi8_mask(a, zero_); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, _mm512_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0 + + tmp_mask = _mm512_cmpgt_epi8_mask(b, zero_); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, _mm512_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(a2, zero_); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, _mm512_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b2, zero_); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, _mm512_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0 + //_mm512_store_si512(&pr[t], d); + _mm512_mask_storeu_epi8(&pr[t], msk, d); + } + else + #endif + { + tmp_mask = _mm512_cmpgt_epi8_mask(a, zero_); + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, _mm512_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b, zero_); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, _mm512_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(a2, zero_); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, _mm512_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b2, zero_); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, _mm512_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0 + // int64_t tic = _rdtsc(); + _mm512_store_si512(&pr[t], d); + // int64_t toc = _rdtsc(); + //fprintf(stderr, "store time: %d %ld\n", t, toc - tic); + } + } +*/ + //continue; + } else { // gap right-alignment + // continue; + __m512i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + // off[r] = stn, off_end[r] = enn; + // fprintf(stderr, "t: %d, st0: %d\n", st_, st0); + #if TIMER_I + uint64_t tic = _rdtsc(); + #endif + for (t = st_; t <= en_; ++t) { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + + //__dp_code_block1_pcl; + __dp_code_block1_pcl_512_mm2_fast; + //d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a),_mm512_set1_epi8(1), zero_); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a), one_, zero_); + z = _mm512_max_epi8(z, a); + // d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b), _mm512_set1_epi8(2), d); // d = z > b? d : 2 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b), two_, d); // d = z > b? d : 2 + z = _mm512_max_epi8(z, b); + // d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a2), _mm512_set1_epi8(3), d); // d = z > a2? d : 3 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a2), three_, d); // d = z > a2? d : 3 + z = _mm512_max_epi8(z, a2); + // d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b2), _mm512_set1_epi8(4), d); // d = z > b2? d : 4 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b2), four_, d); // d = z > b2? d : 4 + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + + // __dp_code_block2_pcl; + + #if 1 + __mmask64 msk; + if (t == en_) { + // __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + off_end[r] -= (4-ind)*16; + + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + else { + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + #endif + + #if 1 + if (t == en_) { + //__mmask64 msk; + //int ind = (en0 - t*64 + 16)/16; + //msk = msk_ar2[ind]; + // off_end[r] -= (4-ind)*16; + + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a, zero_), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, _mm512_set1_epi8(0x08), zero_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b, zero_), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, _mm512_set1_epi8(0x10), zero_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a2); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a2, zero_), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, _mm512_set1_epi8(0x20), zero_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b2); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b2, zero_), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, _mm512_set1_epi8(0x40), zero_)); // d = b > 0? 1<<6 : 0 + // _mm512_store_si512(&pr[t], d); + _mm512_mask_storeu_epi8(&pr[t], msk, d); + + } + else + #endif + { + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a); + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a, zero_), qe_)); + // d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, _mm512_set1_epi8(0x08), zero_)); // d = a > 0? 1<<3 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s1_, zero_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b, zero_), qe_)); + // d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, _mm512_set1_epi8(0x10), zero_)); // d = b > 0? 1<<4 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s2_, zero_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a2); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a2, zero_), qe2_)); + // d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, _mm512_set1_epi8(0x20), zero_)); // d = a > 0? 1<<5 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s3_, zero_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b2); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b2, zero_), qe2_)); + // d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, _mm512_set1_epi8(0x40), zero_)); // d = b > 0? 1<<6 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s4_, zero_)); // d = b > 0? 1<<6 : 0 + _mm512_store_si512(&pr[t], d); + + } + } + #if TIMER_I + uint64_t toc = _rdtsc(); + tprof[B4][0] += toc - tic; + #endif + } +// continue; + + #if TIMER_I + tic = _rdtsc(); + #endif + if (!approx_max) { // find the exact max with a 32-bit score array + int32_t max_H, max_t; + // compute H[], max_H and max_t + if (r > 0) { + int32_t HH[16], tt[16], en1 = st0 + (en0 - st0) / 16 * 16, i; //TODO: saurabh HH[4] -> HH[16]. 4*4?? + __m512i max_H_, max_t_; + max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] : H[en0] + v8[en0]; // special casing the last element + + max_t = en0; + max_H_ = _mm512_set1_epi32(max_H); + max_t_ = _mm512_set1_epi32(max_t); + for (t = st0; t < en1; t += /*4*/16) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t; + __m512i H1, tmp, t_; + __mmask16 tmp_mask; + H1 = _mm512_loadu_si512((__m512i*)&H[t]); + __m128i t__ = _mm_load_si128((__m128i*) &v8[t]); + t_ = _mm512_cvtepi8_epi32(t__); + // t_ = _mm512_setr_epi32(v8[t], v8[t+1], v8[t+2], v8[t+3],v8[t+4], v8[t+5], v8[t+6], v8[t+7],v8[t+8], v8[t+9], v8[t+10], v8[t+11],v8[t+12], v8[t+13], v8[t+14], v8[t+15]); //TODO: need to understand well + H1 = _mm512_add_epi32(H1, t_); + _mm512_storeu_si512((__m512i*)&H[t], H1); + // new code by vasim: making it 4 lanes to match accuracy + { + #if 0 + t_ = _mm512_set1_epi32(t); + __m512i shfH, shft, max1, max2; + for (int l=0; l<4; l++) { + //t_ = _mm512_set1_epi32(t + l*4); + // __m512i shfH = H1, shft = t_; + shfH = _mm512_shuffle_i32x4(H1, H1, l); + shft = _mm512_shuffle_i32x4(t_, t_, l); + tmp_mask = _mm512_cmpgt_epi32_mask(shfH, max_H_); + max_H_ = _mm512_mask_blend_epi32(tmp_mask, max_H_, shfH); + max_t_ = _mm512_mask_blend_epi32(tmp_mask, max_t_, shft); + t_ = _mm512_add_epi32(t_, four32_); + } + #else + __m512i shfH, shft, max1, max2; + t_ = _mm512_set1_epi32(t); + t_ = _mm512_add_epi32(t_, bt32_); + shfH = _mm512_shuffle_i32x4(H1, H1, 0x31); + shft = _mm512_shuffle_i32x4(t_, t_, 0x31); + tmp_mask = _mm512_cmpgt_epi32_mask(shfH, H1); + max1 = _mm512_mask_blend_epi32(tmp_mask, H1, shfH); + max2 = _mm512_mask_blend_epi32(tmp_mask, t_, shft); + shfH = _mm512_shuffle_i32x4(max1, max1, 0x2); + shft = _mm512_shuffle_i32x4(max2, max2, 0x2); + tmp_mask = _mm512_cmpgt_epi32_mask(shfH, max1); + max1 = _mm512_mask_blend_epi32(tmp_mask, max1, shfH); + max2 = _mm512_mask_blend_epi32(tmp_mask, max2, shft); + tmp_mask = _mm512_cmpgt_epi32_mask(max1, max_H_); + max_H_ = _mm512_mask_blend_epi32(tmp_mask, max_H_, max1); + max_t_ = _mm512_mask_blend_epi32(tmp_mask, max_t_, max2); + #endif + } + } + _mm512_storeu_si512((__m512i*)HH, max_H_); + _mm512_storeu_si512((__m512i*)tt, max_t_); + + // fprintf(stderr, "max_H: %d, max_t: %d\n", max_H, max_t); + int rem = (en0 - t) / 4; + for (int l=0; l HH[j]) { + HH[j] = H[t]; + tt[j] = bt; + } + t++; + } + } + #pragma unroll + for (i = 0; i < 4; ++i) + if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i; + + for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE + H[t] += (int32_t)v8[t]; + if (H[t] > max_H) { + max_H = H[t], max_t = t; + } + } + } else H[0] = v8[0] - qe, max_H = H[0], max_t = 0; // special casing r==0 + + // update ez + if (en0 == tlen - 1 && H[en0] > ez->mte) { + // ez->mte = H[en0], ez->mte_q = r - en; + ez->mte = H[en0], ez->mte_q = r - en_new; + // fprintf(stderr, "r: %d, en: %d\n", r, en); + } + if (r - st0 == qlen - 1 && H[st0] > ez->mqe) { + ez->mqe = H[st0], ez->mqe_t = st0; + } + + if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e2)) { + // fprintf(stderr, "Dropping out, r: %d\n", r); + break; + } + if (r == qlen + tlen - 2 && en0 == tlen - 1) { + ez->score = H[tlen - 1]; + // fprintf(stderr, "ez->score: %d\n", ez->score); + } + } else { // find approximate max; Z-drop might be inaccurate, too. + if (r > 0) { + if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) { + int32_t d0 = v8[last_H0_t]; + int32_t d1 = u8[last_H0_t + 1]; + if (d0 > d1) H0 += d0; + else H0 += d1, ++last_H0_t; + // fprintf(stderr, "1. r: %d, d0: %d, d1: %d H0: %d, last_H0_t: %d\n", r, d0, d1, H0, last_H0_t); + } else if (last_H0_t >= st0 && last_H0_t <= en0) { + H0 += v8[last_H0_t]; + // fprintf(stderr, "2. r: %d, H0: %d, last_H0_t: %d, v8: %d\n", r, H0, last_H0_t, v8[last_H0_t]); + } else { + ++last_H0_t, H0 += u8[last_H0_t]; + // fprintf(stderr, "3. r: %d, H0: %d, last_H0_t: %d, u8: %d\n", r, H0, last_H0_t, u8[last_H0_t]); + } + } else H0 = v8[0] - qe, last_H0_t = 0; + if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e2)) { + // fprintf(stderr, "Dropping II\n"); + break; + } + if (r == qlen + tlen - 2 && en0 == tlen - 1) { + ez->score = H0; + // fprintf(stderr, ">>r: %d, ez->score: %d\n", r, ez->score); + } + } + // last_st = st, last_en = en; + last_st = st_new, last_en = en_new; + + #if TIMER_I + toc = _rdtsc(); + tprof[B5][0] += toc - tic; + #endif + } + #if TIMER_O + uint64_t toc = _rdtsc(); + tprof[OUTER][0] += toc - tic; + #endif + + kfree(km, mem); + if (!approx_max) kfree(km, H); + + #if 1 + if (with_cigar) { // backtrack + int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR); + if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) { + ez->reach_end = 1; + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (ez->max_t >= 0 && ez->max_q >= 0) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } + kfree(km, mem2); kfree(km, off); + } + #endif + + //kfree_all(); +} + + + +void ksw_extd2_avx512_v1(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) +{ + __m512i bt32_ = _mm512_setr_epi32(0,0,0,0,4,4,4,4,8,8,8,8,12,12,12,12); + + int8_t index[64] __attribute((aligned(64))); + + for (int i=0; i<64; i++) + index[i] = i%16 - 1; + + index[0] = 15; + index[16] = 31; + index[32] = 47; + index[48] = 63; + + __m512i shf512a, shf512b, slli512; + __m512i ind512_slli = _mm512_load_si512((__m512i*) index); + __mmask8 mska = 0x90; + __mmask64 mskb = 0x0001000100010000; + __mmask64 mskc = 0x1; + __mmask64 mskc_ar[4] = {0x1, 0x10000, 0x100000000, 0x1000000000000}; + + + + #define __dp_code_block1_pcl \ + /*__mmask64 mskc_ = mskc; */ \ + /*if (t == st_) */ \ + __mmask64 mskc_ = (t == st_) ? mskc_ar[(st0 - t*64)/16]:mskc; \ + z = _mm512_load_si512(&s[t]); \ + xt1 = _mm512_load_si512(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \ + /* tmp = _mm_srli_si128(xt1, 15); */ /* tmp <- x[r-1][t+15] */ \ + tmp = _mm512_set1_epi8(((int8_t*)x)[t*64 + 63]); \ + /* xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_);*/ /* xt1 <- x[r-1][t-1..t+14] */ \ + shf512a = _mm512_shuffle_epi8(xt1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + xt1 = _mm512_mask_blend_epi8(mskc_, slli512, x1_); \ + x1_ = tmp; \ + vt1 = _mm512_load_si512(&v[t]); /* vt1 <- v[r-1][t..t+15] */ \ + /* tmp = _mm_srli_si128(vt1, 15); */ /* tmp <- v[r-1][t+15] */ \ + tmp = _mm512_set1_epi8(((int8_t*)v)[t*64 + 63]); \ + /* vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); *//* vt1 <- v[r-1][t-1..t+14] */ \ + shf512a = _mm512_shuffle_epi8(vt1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + vt1 = _mm512_mask_blend_epi8(mskc_, slli512, v1_); \ + v1_ = tmp; \ + a = _mm512_add_epi8(xt1, vt1); /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \ + ut = _mm512_load_si512(&u[t]); /* ut <- u[t..t+15] */ \ + b = _mm512_add_epi8(_mm512_load_si512(&y[t]), ut); /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */ \ + x2t1= _mm512_load_si512(&x2[t]); \ + /* tmp = _mm_srli_si128(x2t1, 15);*/ \ + tmp = _mm512_set1_epi8(((int8_t*)x2)[t*64 + 63]); \ + /* x2t1= _mm_or_si128(_mm_slli_si128(x2t1, 1), x21_); */ \ + shf512a = _mm512_shuffle_epi8(x2t1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + x2t1 = _mm512_mask_blend_epi8(mskc_, slli512, x21_); \ + x21_= tmp; \ + a2= _mm512_add_epi8(x2t1, vt1); \ + b2= _mm512_add_epi8(_mm512_load_si512(&y2[t]), ut); + + + #define __dp_code_block2_pcl \ + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \ + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \ + tmp = _mm512_sub_epi8(z, q_); \ + a = _mm512_sub_epi8(a, tmp); \ + b = _mm512_sub_epi8(b, tmp); \ + tmp = _mm512_sub_epi8(z, q2_); \ + a2= _mm512_sub_epi8(a2, tmp); \ + b2= _mm512_sub_epi8(b2, tmp); + + __mmask64 msk_ar2[5] = {0xFFFF, 0xFFFF, 0xFFFFFFFF, 0xFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF}; + + int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc, long_thres, long_diff; + int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX); + int32_t *H = 0, H0 = 0, last_H0_t = 0; + uint8_t *qr, *sf, *mem, *mem2 = 0; + + __m512i q_, q2_, qe_, qe2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_; + __m512i *u, *v, *x, *y, *x2, *y2, *s, *p = 0; + __m512i one_, two_, three_, four_, s1_, s2_, s3_, s4_; + + ksw_reset_extz(ez); + if (m <= 1 || qlen <= 0 || tlen <= 0) return; + + if (q2 + e2 < q + e) t = q, q = q2, q2 = t, t = e, e = e2, e2 = t; // make sure q+e no larger than q2+e2 + s1_ = _mm512_set1_epi8(0x08); + s2_ = _mm512_set1_epi8(0x10); + s3_ = _mm512_set1_epi8(0x20); + s4_ = _mm512_set1_epi8(0x40); + + one_ = _mm512_set1_epi8(1); + two_ = _mm512_set1_epi8(2); + three_ = _mm512_set1_epi8(3); + four_ = _mm512_set1_epi8(4); + + zero_ = _mm512_set1_epi8(0); + q_ = _mm512_set1_epi8(q); + q2_ = _mm512_set1_epi8(q2); + qe_ = _mm512_set1_epi8(q + e); + qe2_ = _mm512_set1_epi8(q2 + e2); + sc_mch_ = _mm512_set1_epi8(mat[0]); + sc_mis_ = _mm512_set1_epi8(mat[1]); + sc_N_ = mat[m*m-1] == 0? _mm512_set1_epi8(-e2) : _mm512_set1_epi8(mat[m*m-1]); + m1_ = _mm512_set1_epi8(m - 1); // wildcard + + if (w < 0) w = tlen > qlen? tlen : qlen; + wl = wr = w; + tlen_ = (tlen + 63) / 64; + n_col_ = qlen < tlen? qlen : tlen; + n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 63) / 64 + 1; + qlen_ = (qlen + 63) / 64; + for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) { + max_sc = max_sc > mat[t]? max_sc : mat[t]; + min_sc = min_sc < mat[t]? min_sc : mat[t]; + } + if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches + + long_thres = e != e2? (q2 - q) / (e - e2) - 1 : 0; + if (q2 + e2 + long_thres * e2 > q + e + long_thres * e) + ++long_thres; + long_diff = long_thres * (e - e2) - (q2 - q) - e2; + + mem = (uint8_t*)kcalloc(km, tlen_ * 8 + qlen_ + 1 + 63, 64); + u = (__m512i*)(((size_t)mem + 63) >> 6 << 6); // 16-byte aligned + v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_; + s = y2 + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 64; + memset(u, -q - e, tlen_ * 64); + memset(v, -q - e, tlen_ * 64); + memset(x, -q - e, tlen_ * 64); + memset(y, -q - e, tlen_ * 64); + memset(x2, -q2 - e2, tlen_ * 64); + memset(y2, -q2 - e2, tlen_ * 64); + if (!approx_max) { + H = (int32_t*)kmalloc(km, tlen_ * 64 * 4); + for (t = 0; t < tlen_ * 64; ++t) H[t] = KSW_NEG_INF; + } + if (with_cigar) { + mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 64); + p = (__m512i*)(((size_t)mem2 + 63) >> 6 << 6); + off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2); + off_end = off + qlen + tlen - 1; + } + + for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t]; + memcpy(sf, target, tlen); + + + for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) { + int st = 0, en = tlen - 1, st0, en0, st_, en_; + int8_t x1, x21, v1; + uint8_t *qrr = qr + (qlen - 1 - r); + int8_t *u8 = (int8_t*)u, *v8 = (int8_t*)v, *x8 = (int8_t*)x, *x28 = (int8_t*)x2; + __m512i x1_, x21_, v1_; + // find the boundaries + if (st < r - qlen + 1) st = r - qlen + 1; + if (en > r) en = r; + if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil + if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor + if (st > en) { + ez->zdropped = 1; + break; + } + st0 = st, en0 = en; + int st_new = st / 16 * 16, en_new = (en + 16) / 16 * 16 - 1; + // int st_new = st / 64 * 64, en_new = (en + 64) / 64 * 64 - 1; + //int stb = st, enb = en; + st = st / 64 * 64, en = (en + 64) / 64 * 64 - 1; + //int stn = stb / 16 * 16, enn = (enb + 16) / 16 * 16 - 1; + // set boundary conditions + if (st_new > 0) { + if (st_new - 1 >= last_st && st_new - 1 <= last_en) { + x1 = x8[st_new - 1], x21 = x28[st_new - 1], v1 = v8[st_new - 1]; // (r-1,s-1) calculated in the last round + } else { + x1 = -q - e, x21 = -q2 - e2; + v1 = -q - e; + } + } else { + x1 = -q - e, x21 = -q2 - e2; + v1 = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; + } + if (en_new >= r) { + ((int8_t*)y)[r] = -q - e, ((int8_t*)y2)[r] = -q2 - e2; + u8[r] = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; + } + + // loop fission: set scores first + if (!(flag & KSW_EZ_GENERIC_SC)) { + for (t = st0; t <= en0; t += 64) { + __m512i sq, st, tmp_512; + __mmask64 tmp, mask; + sq = _mm512_loadu_si512((__m512i*)&sf[t]); + st = _mm512_loadu_si512((__m512i*)&qrr[t]); +// mask = _mm512_or_si512(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_)); + mask = (_mm512_cmpeq_epi8_mask(sq, m1_) | _mm512_cmpeq_epi8_mask(st, m1_)); + tmp = _mm512_cmpeq_epi8_mask(sq, st); + + tmp_512 = _mm512_mask_blend_epi8(tmp, sc_mis_, sc_mch_); + tmp_512 = _mm512_mask_blend_epi8(mask, tmp_512, sc_N_); + if (t + 64 > en0) + { + __mmask64 msk; + int ind = (en0 - t + 16)/16; + //assert(ind >= 0 && ind < 5); + msk = msk_ar2[ind]; + _mm512_mask_storeu_epi8((__m512i*)((int8_t*)s + t), msk, tmp_512); + } + else + _mm512_storeu_si512((__m512i*)((int8_t*)s + t), tmp_512); + + } + } else { + for (t = st0; t <= en0; ++t) + ((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]]; + } + + // core loop + // fprintf(stderr, "- r: %d, x1: %d, x21: %d, v1: %d, en_new: %d, e: %d, q: %d\n", + //r, x1, x21, v1, en_new, e, q); + x1_ = _mm512_set1_epi8((uint8_t)x1); + x21_ = _mm512_set1_epi8((uint8_t)x21); + v1_ = _mm512_set1_epi8((uint8_t)v1); + + //st_ = st / 16, en_ = en / 16; + st_ = st / 64, en_ = en / 64; + //assert(en_ - st_ + 1 <= n_col_); + if (!with_cigar) { // score only + for (t = st_; t <= en_; ++t) { + __m512i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __dp_code_block1_pcl; + + z = _mm512_max_epi8(z, a); + z = _mm512_max_epi8(z, b); + z = _mm512_max_epi8(z, a2); + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; // save u[] and v[]; update a, b, a2 and b2 + if (t == en_) { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + // fprintf(stderr, "en0: %d, t: %d, ind: %d, msk: %d\n", en0, t, ind, msk); + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + + } + else { + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + + if (t == en_) { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + //assert(ind >= 0); + msk = msk_ar2[ind]; + // fprintf(stderr, "en0: %d, t: %d, ind: %d, msk: %d\n", en0, t, ind, msk); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_max_epi8(a, zero_), qe_)); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_max_epi8(b, zero_), qe_)); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_max_epi8(a2, zero_), qe2_)); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_max_epi8(b2, zero_), qe2_)); + } + else + { + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_max_epi8(a, zero_), qe_)); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_max_epi8(b, zero_), qe_)); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_max_epi8(a2, zero_), qe2_)); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_max_epi8(b2, zero_), qe2_)); + } + // for (int l=0; l<64; l++) + // fprintf(stderr, "%d ", ((int8_t*)x)[l]); + + } + + } else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment + __m512i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + + + for (t = st_; t < en_; ++t) { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a, z),zero_, one_); // d = a > z? 1 : 0 + z = _mm512_max_epi8(z, a); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b, z), d, two_); // d = b > z? 2 : d + z = _mm512_max_epi8(z, b); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a2, z), d, three_); // d = a2 > z? 3 : d + z = _mm512_max_epi8(z, a2); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b2, z), d, four_); // d = b2 > z? 4 : d + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + tmp_mask = _mm512_cmpgt_epi8_mask(a, zero_); + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s1_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b, zero_); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s2_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(a2, zero_); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s3_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b2, zero_); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s4_)); // d = b > 0? 1<<6 : 0 + _mm512_store_si512(&pr[t], d); + } + { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a, z),zero_, one_); // d = a > z? 1 : 0 + z = _mm512_max_epi8(z, a); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b, z), d, two_); // d = b > z? 2 : d + z = _mm512_max_epi8(z, b); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a2, z), d, three_); // d = a2 > z? 3 : d + z = _mm512_max_epi8(z, a2); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b2, z), d, four_); // d = b2 > z? 3 : d + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + // //assert(ind >= 0 && ind < 5); + msk = msk_ar2[ind]; + + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + + { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + off_end[r] -= (4-ind)*16; + + tmp_mask = _mm512_cmpgt_epi8_mask(a, zero_); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s1_)); // d = a > 0? 1<<3 : 0 + + tmp_mask = _mm512_cmpgt_epi8_mask(b, zero_); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s2_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(a2, zero_); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s3_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b2, zero_); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s4_)); // d = b > 0? 1<<6 : 0 + //_mm512_store_si512(&pr[t], d); + _mm512_mask_storeu_epi8(&pr[t], msk, d); + } + } + + + } else { // gap right-alignment + __m512i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + // off[r] = stn, off_end[r] = enn; + // fprintf(stderr, "t: %d, st0: %d\n", st_, st0); + for (t = st_; t < en_; ++t) { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a), one_, zero_); + z = _mm512_max_epi8(z, a); + // d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b), _mm512_set1_epi8(2), d); + // d = z > b? d : 2 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b), two_, d); // d = z > b? d : 2 + z = _mm512_max_epi8(z, b); + // d = z > a2? d : 3 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a2), three_, d); // d = z > a2? d : 3 + z = _mm512_max_epi8(z, a2); + // d = z > b2? d : 4 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b2), four_, d); // d = z > b2? d : 4 + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + + // __dp_code_block2_pcl; + + { + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + + { + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a); + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a, zero_), qe_)); + // d = a > 0? 1<<3 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s1_, zero_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b, zero_), qe_)); + // d = b > 0? 1<<4 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s2_, zero_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a2); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a2, zero_), qe2_)); + // d = a > 0? 1<<5 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s3_, zero_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b2); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b2, zero_), qe2_)); + // d = b > 0? 1<<6 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s4_, zero_)); // d = b > 0? 1<<6 : 0 + _mm512_store_si512(&pr[t], d); + + } + } + //for (t = st_; t <= en_; ++t)// Last iteration unrolled + { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a), one_, zero_); + z = _mm512_max_epi8(z, a); + // d = z > b? d : 2 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b), two_, d); // d = z > b? d : 2 + z = _mm512_max_epi8(z, b); + // d = z > a2? d : 3 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a2), three_, d); // d = z > a2? d : 3 + z = _mm512_max_epi8(z, a2); + // d = z > b2? d : 4 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b2), four_, d); // d = z > b2? d : 4 + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + + // __dp_code_block2_pcl; + + __mmask64 msk; + { + // __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + off_end[r] -= (4-ind)*16; + + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + + { + + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a, zero_), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s1_, zero_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b, zero_), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s2_, zero_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a2); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a2, zero_), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s3_, zero_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b2); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b2, zero_), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s4_, zero_)); // d = b > 0? 1<<6 : 0 + // _mm512_store_si512(&pr[t], d); + _mm512_mask_storeu_epi8(&pr[t], msk, d); + + } + } + } + + if (!approx_max) { // find the exact max with a 32-bit score array + int32_t max_H, max_t; + // compute H[], max_H and max_t + if (r > 0) { + int32_t HH[16], tt[16], en1 = st0 + (en0 - st0) / 16 * 16, i; + __m512i max_H_, max_t_; + max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] : H[en0] + v8[en0]; // special casing the last element + + max_t = en0; + max_H_ = _mm512_set1_epi32(max_H); + max_t_ = _mm512_set1_epi32(max_t); + for (t = st0; t < en1; t += /*4*/16) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t; + __m512i H1, t_; + __mmask16 tmp_mask; + H1 = _mm512_loadu_si512((__m512i*)&H[t]); + __m128i t__ = _mm_load_si128((__m128i*) &v8[t]); + t_ = _mm512_cvtepi8_epi32(t__); + H1 = _mm512_add_epi32(H1, t_); + _mm512_storeu_si512((__m512i*)&H[t], H1); + // making it 4 lanes to match accuracy + + __m512i shfH, shft, max1, max2; + t_ = _mm512_set1_epi32(t); + t_ = _mm512_add_epi32(t_, bt32_); + shfH = _mm512_shuffle_i32x4(H1, H1, 0x31); + shft = _mm512_shuffle_i32x4(t_, t_, 0x31); + tmp_mask = _mm512_cmpgt_epi32_mask(shfH, H1); + max1 = _mm512_mask_blend_epi32(tmp_mask, H1, shfH); + max2 = _mm512_mask_blend_epi32(tmp_mask, t_, shft); + shfH = _mm512_shuffle_i32x4(max1, max1, 0x2); + shft = _mm512_shuffle_i32x4(max2, max2, 0x2); + tmp_mask = _mm512_cmpgt_epi32_mask(shfH, max1); + max1 = _mm512_mask_blend_epi32(tmp_mask, max1, shfH); + max2 = _mm512_mask_blend_epi32(tmp_mask, max2, shft); + tmp_mask = _mm512_cmpgt_epi32_mask(max1, max_H_); + max_H_ = _mm512_mask_blend_epi32(tmp_mask, max_H_, max1); + max_t_ = _mm512_mask_blend_epi32(tmp_mask, max_t_, max2); + + + } + _mm512_storeu_si512((__m512i*)HH, max_H_); + _mm512_storeu_si512((__m512i*)tt, max_t_); + + int rem = (en0 - t) / 4; + for (int l=0; l HH[j]) { + HH[j] = H[t]; + tt[j] = bt; + } + t++; + } + } + + for (i = 0; i < 4; ++i) + if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i; + + for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE + H[t] += (int32_t)v8[t]; + if (H[t] > max_H) { + max_H = H[t], max_t = t; + } + } + } else H[0] = v8[0] - qe, max_H = H[0], max_t = 0; // special casing r==0 + + // update ez + if (en0 == tlen - 1 && H[en0] > ez->mte) { + ez->mte = H[en0], ez->mte_q = r - en_new; + } + if (r - st0 == qlen - 1 && H[st0] > ez->mqe) { + ez->mqe = H[st0], ez->mqe_t = st0; + } + + if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e2)) { + break; + } + if (r == qlen + tlen - 2 && en0 == tlen - 1) { + ez->score = H[tlen - 1]; + } + } else { // find approximate max; Z-drop might be inaccurate, too. + if (r > 0) { + if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) { + int32_t d0 = v8[last_H0_t]; + int32_t d1 = u8[last_H0_t + 1]; + if (d0 > d1) H0 += d0; + else H0 += d1, ++last_H0_t; + } else if (last_H0_t >= st0 && last_H0_t <= en0) { + H0 += v8[last_H0_t]; + } else { + ++last_H0_t, H0 += u8[last_H0_t]; + } + } else H0 = v8[0] - qe, last_H0_t = 0; + if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e2)) { + break; + } + if (r == qlen + tlen - 2 && en0 == tlen - 1) { + ez->score = H0; + } + } + // last_st = st, last_en = en; + last_st = st_new, last_en = en_new; + + } + + kfree(km, mem); + if (!approx_max) kfree(km, H); + if (with_cigar) { // backtrack + int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR); + if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) { + ez->reach_end = 1; + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (ez->max_t >= 0 && ez->max_q >= 0) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } + kfree(km, mem2); kfree(km, off); + } +#undef __dp_code_block1_pcl +#undef __dp_code_block2_pcl +} + +#endif + +#ifdef __AVX2__ +/* +__m256i get_mask_store(__m256i msk, void* addr){ + return _mm256_or_si256(msk , _mm256_loadu_si256((__m256i*)addr)); +} +*/ +__m256i get_mask_store(__m256i msk, void* addr, __m256i store_data){ + return _mm256_blendv_epi8(_mm256_loadu_si256((__m256i*)addr), store_data, msk); +} + +extern void *km1; +extern uint64_t km_size;// = 500000000; // 500 MB +extern int km_top; +extern uint64_t avg; +extern uint64_t alignment_time; +void *kcalloc_(void* km, int count, int size) +{ + +// assert(km_top == 0); +// assert(count*size < km_size); + if(!(km_top == 0 && count*size < km_size)){ + fprintf(stderr, "%d %lld %lld %lld %lld", km_top, count, size, km_size, count*size); + + } + assert(km_top == 0); + assert(count*size < km_size); + km_top += count*size + 1024; + memset(km, 0, count * size); + + // printf("km_top: %d\n", km_top); + return km; +} + +void *kmalloc_(void* km, int count) { + if(km_top + count >= km_size) + printf("count: %d\n", count); + assert(km_top + count < km_size); + void *mem = (void*) ((int8_t*) km + km_top); + km_top += count + 1024; + // printf("km_top: %d\n", km_top); + return mem; +} + +void kfree_all() { km_top = 0;} + +void ksw_extd2_avx2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) +{ + + + __m256i bt32_ = _mm256_setr_epi32(0,0,0,0,4,4,4,4);//8,8,8,8,12,12,12,12); + + int8_t index[32] __attribute((aligned(64))); + + for (int i=0; i<32; i++) + index[i] = i%16 - 1; + + index[0] = 15; + index[16] = 31; +// index[32] = 47; +// index[48] = 63; + + __m256i shf256a, shf256b, slli256; + __m256i ind256_slli = _mm256_load_si256((__m256i*) index); + //__mmask8 mska = 0x00;//0x90 + //__mmask32 mskb = 0x00010000;//0000 0000 0000 0001 0000 0000 0000 0000 + //__mmask32 mskc = 0x1; + + __m256i mskb_v = _mm256_set_epi32(0,0,0,255,0,0,0,0); + __m256i mskc_v = _mm256_set_epi32(0,0,0,0,0,0,0,255); + + //__mmask32 mskc_ar[2] = {0x1, 0x10000}; + __m256i mskc_ar_v[2];// = {0x1, 0x10000}; + mskc_ar_v[0] = _mm256_set_epi32(0,0,0,0,0,0,0,255); + mskc_ar_v[1] = _mm256_set_epi32(0,0,0,255,0,0,0,0); + + + + #define __dp_code_block1_pcl \ + /* __mmask32 mskc_ = (t == st_) ? mskc_ar[(st0 - t*32)/16]:mskc; */ \ + __m256i mskc_ = (t == st_) ? mskc_ar_v[(st0 - t*32)/16]:mskc_v; \ + z = _mm256_load_si256(&s[t]); \ + xt1 = _mm256_load_si256(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \ + /* tmp = _mm_srli_si128(xt1, 15); */ /* tmp <- x[r-1][t+15] */ \ + tmp = _mm256_set1_epi8(((int8_t*)x)[t*32 + 31]); \ + /* xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_);*/ /* xt1 <- x[r-1][t-1..t+14] */ \ + shf256a = _mm256_shuffle_epi8(xt1, ind256_slli); \ + /*shf256b = _mm256_shuffle_i32x4(shf256a, shf256a, 0x00); */ \ + shf256b = _mm256_permute2x128_si256(shf256a, shf256a, 0); \ +/* slli256 = _mm256_mask_blend_epi8(mskb, shf256a, shf256b); */ \ + slli256 = _mm256_blendv_epi8(shf256a, shf256b, mskb_v); \ + xt1 = _mm256_blendv_epi8(slli256, x1_, mskc_); \ + x1_ = tmp; \ + vt1 = _mm256_load_si256(&v[t]); /* vt1 <- v[r-1][t..t+15] */ \ + /* tmp = _mm_srli_si128(vt1, 15); */ /* tmp <- v[r-1][t+15] */ \ + tmp = _mm256_set1_epi8(((int8_t*)v)[t*32 + 31]); \ + /* vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); *//* vt1 <- v[r-1][t-1..t+14] */ \ + shf256a = _mm256_shuffle_epi8(vt1, ind256_slli); \ + /*shf256b = _mm256_shuffle_i32x4(shf256a, shf256a, 0x00); */ \ + shf256b = _mm256_permute2x128_si256(shf256a, shf256a, 0); \ +/* slli256 = _mm256_mask_blend_epi8(mskb, shf256a, shf256b); */ \ + slli256 = _mm256_blendv_epi8(shf256a, shf256b, mskb_v); \ + vt1 = _mm256_blendv_epi8(slli256, v1_, mskc_); \ + v1_ = tmp; \ + a = _mm256_add_epi8(xt1, vt1); /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \ + ut = _mm256_load_si256(&u[t]); /* ut <- u[t..t+15] */ \ + b = _mm256_add_epi8(_mm256_load_si256(&y[t]), ut); /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */ \ + x2t1= _mm256_load_si256(&x2[t]); \ + /* tmp = _mm_srli_si128(x2t1, 15);*/ \ + tmp = _mm256_set1_epi8(((int8_t*)x2)[t*32 + 31]); \ + /* x2t1= _mm_or_si128(_mm_slli_si128(x2t1, 1), x21_); */ \ + shf256a = _mm256_shuffle_epi8(x2t1, ind256_slli); \ + shf256b = _mm256_permute2x128_si256(shf256a, shf256a, 0); \ + slli256 = _mm256_blendv_epi8(shf256a, shf256b, mskb_v); \ + x2t1 = _mm256_blendv_epi8(slli256, x21_, mskc_); \ + x21_= tmp; \ + a2= _mm256_add_epi8(x2t1, vt1); \ + b2= _mm256_add_epi8(_mm256_load_si256(&y2[t]), ut); + + + #define __dp_code_block2_pcl \ + _mm256_storeu_si256(&u[t], _mm256_sub_epi8(z, vt1)); /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \ + _mm256_storeu_si256(&v[t], _mm256_sub_epi8(z, ut)); /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \ + tmp = _mm256_sub_epi8(z, q_); \ + a = _mm256_sub_epi8(a, tmp); \ + b = _mm256_sub_epi8(b, tmp); \ + tmp = _mm256_sub_epi8(z, q2_); \ + a2= _mm256_sub_epi8(a2, tmp); \ + b2= _mm256_sub_epi8(b2, tmp); + + //__mmask32 msk_ar2[3] = {0xFFFF, 0xFFFF, 0xFFFFFFFF}; + __m256i msk_ar2_v[3]; + msk_ar2_v[0] = _mm256_set_epi32(0,0,0,0,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF); + msk_ar2_v[1] = _mm256_set_epi32(0,0,0,0,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF); + msk_ar2_v[2] = _mm256_set_epi32(0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF,0xFFFFFFFF); + + + + int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc, long_thres, long_diff; + int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX); + int32_t *H = 0, H0 = 0, last_H0_t = 0; + uint8_t *qr, *sf, *mem, *mem2 = 0; + + __m256i q_, q2_, qe_, qe2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_; + __m256i *u, *v, *x, *y, *x2, *y2, *s, *p = 0; + __m256i one_, two_, three_, four_, s1_, s2_, s3_, s4_; + + ksw_reset_extz(ez); + if (m <= 1 || qlen <= 0 || tlen <= 0) return; + + if (q2 + e2 < q + e) t = q, q = q2, q2 = t, t = e, e = e2, e2 = t; // make sure q+e no larger than q2+e2 + s1_ = _mm256_set1_epi8(0x08); + s2_ = _mm256_set1_epi8(0x10); + s3_ = _mm256_set1_epi8(0x20); + s4_ = _mm256_set1_epi8(0x40); + + one_ = _mm256_set1_epi8(1); + two_ = _mm256_set1_epi8(2); + three_ = _mm256_set1_epi8(3); + four_ = _mm256_set1_epi8(4); + + zero_ = _mm256_set1_epi8(0); + q_ = _mm256_set1_epi8(q); + q2_ = _mm256_set1_epi8(q2); + qe_ = _mm256_set1_epi8(q + e); + qe2_ = _mm256_set1_epi8(q2 + e2); + sc_mch_ = _mm256_set1_epi8(mat[0]); + sc_mis_ = _mm256_set1_epi8(mat[1]); + sc_N_ = mat[m*m-1] == 0? _mm256_set1_epi8(-e2) : _mm256_set1_epi8(mat[m*m-1]); + m1_ = _mm256_set1_epi8(m - 1); // wildcard + + if (w < 0) w = tlen > qlen? tlen : qlen; + wl = wr = w; + tlen_ = (tlen + 31) / 32; + n_col_ = qlen < tlen? qlen : tlen; + n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 31) / 32 + 1; + qlen_ = (qlen + 31) / 32; + for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) { + max_sc = max_sc > mat[t]? max_sc : mat[t]; + min_sc = min_sc < mat[t]? min_sc : mat[t]; + } + if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches + + long_thres = e != e2? (q2 - q) / (e - e2) - 1 : 0; + if (q2 + e2 + long_thres * e2 > q + e + long_thres * e) + ++long_thres; + long_diff = long_thres * (e - e2) - (q2 - q) - e2; + + //mem = (uint8_t*)kcalloc_(km1, tlen_ * 8 + qlen_ + 1 + 63, 64); + mem = (uint8_t*)kcalloc(km, tlen_ * 8 + qlen_ + 1 + 63, 64); + u = (__m256i*)(((size_t)mem + 31) >> 5 << 5); // 16-byte aligned + v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_; + s = y2 + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 32; + memset(u, -q - e, tlen_ * 32); + memset(v, -q - e, tlen_ * 32); + memset(x, -q - e, tlen_ * 32); + memset(y, -q - e, tlen_ * 32); + memset(x2, -q2 - e2, tlen_ * 32); + memset(y2, -q2 - e2, tlen_ * 32); + if (!approx_max) { + //H = (int32_t*)kmalloc_(km1, tlen_ * 32 * 4); + H = (int32_t*)kmalloc(km, tlen_ * 32 * 4); + for (t = 0; t < tlen_ * 32; ++t) H[t] = KSW_NEG_INF; + } + if (with_cigar) { + //mem2 = (uint8_t*)kmalloc_(km1, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 32); + mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 32); + p = (__m256i*)(((size_t)mem2 + 31) >> 5 << 5); + //off = (int*)kmalloc_(km1, (qlen + tlen - 1) * sizeof(int) * 2); + off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2); + off_end = off + qlen + tlen - 1; + } + for(uint64_t itr = 0; itr < ((qlen + tlen - 1) * n_col_ + 1) * 1; itr ++){ + avg+= mem2[itr]; + } +#ifdef MANUAL_PROFILING + uint64_t align_start = __rdtsc(); +#endif + + + for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t]; + memcpy(sf, target, tlen); + + + for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) { + int st = 0, en = tlen - 1, st0, en0, st_, en_; + int8_t x1, x21, v1; + uint8_t *qrr = qr + (qlen - 1 - r); + int8_t *u8 = (int8_t*)u, *v8 = (int8_t*)v, *x8 = (int8_t*)x, *x28 = (int8_t*)x2; + __m256i x1_, x21_, v1_; + // find the boundaries + if (st < r - qlen + 1) st = r - qlen + 1; + if (en > r) en = r; + if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil + if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor + if (st > en) { + ez->zdropped = 1; + break; + } + st0 = st, en0 = en; + int st_new = st / 16 * 16, en_new = (en + 16) / 16 * 16 - 1; + // int st_new = st / 64 * 64, en_new = (en + 64) / 64 * 64 - 1; + //int stb = st, enb = en; + st = st / 32 * 32, en = (en + 32) / 32 * 32 - 1; + //int stn = stb / 16 * 16, enn = (enb + 16) / 16 * 16 - 1; + // set boundary conditions + if (st_new > 0) { + if (st_new - 1 >= last_st && st_new - 1 <= last_en) { + x1 = x8[st_new - 1], x21 = x28[st_new - 1], v1 = v8[st_new - 1]; // (r-1,s-1) calculated in the last round + } else { + x1 = -q - e, x21 = -q2 - e2; + v1 = -q - e; + } + } else { + x1 = -q - e, x21 = -q2 - e2; + v1 = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; + } + if (en_new >= r) { + ((int8_t*)y)[r] = -q - e, ((int8_t*)y2)[r] = -q2 - e2; + u8[r] = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; + } + + // loop fission: set scores first + if (!(flag & KSW_EZ_GENERIC_SC)) { + for (t = st0; t <= en0; t += 32) { + __m256i sq, st, tmp_256, mask_256; + sq = _mm256_loadu_si256((__m256i*)&sf[t]); + st = _mm256_loadu_si256((__m256i*)&qrr[t]); +// mask = (_mm256_cmpeq_epi8_mask(sq, m1_) | _mm256_cmpeq_epi8_mask(st, m1_)); + mask_256 = _mm256_or_si256(_mm256_cmpeq_epi8(sq, m1_), _mm256_cmpeq_epi8(st, m1_)); + + tmp_256 = _mm256_cmpeq_epi8(sq, st); + + tmp_256 = _mm256_blendv_epi8(sc_mis_, sc_mch_, tmp_256); + tmp_256 = _mm256_blendv_epi8(tmp_256, sc_N_, mask_256); + if (t + 32 > en0) + { + int ind = (en0 - t + 16)/16; + //assert(ind >= 0 && ind < 3); + __m256i msk_v = msk_ar2_v[ind]; + //__m256i str = (get_mask_store(msk_v, ((int8_t*)s + t), tmp_256));// msk_ar2_v[ind]; + //__m256i str =_mm256_and_si256(get_mask_store(msk_v,(int8_t*)s + t), tmp_256);// msk_ar2_v[ind]; + + + //_mm256_storeu_si256((__m256i*)((int8_t*)s + t), str); + + _mm256_maskstore_epi32((int32_t *)((int8_t*)s + t), msk_v, tmp_256); + } + else + _mm256_storeu_si256((__m256i*)((int8_t*)s + t), tmp_256); + + } + } else { + for (t = st0; t <= en0; ++t) + ((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]]; + } + + // core loop + // fprintf(stderr, "- r: %d, x1: %d, x21: %d, v1: %d, en_new: %d, e: %d, q: %d\n", + //r, x1, x21, v1, en_new, e, q); + x1_ = _mm256_set1_epi8((uint8_t)x1); + x21_ = _mm256_set1_epi8((uint8_t)x21); + v1_ = _mm256_set1_epi8((uint8_t)v1); + + //st_ = st / 16, en_ = en / 16; + st_ = st / 32, en_ = en / 32; + //assert(en_ - st_ + 1 <= n_col_); + if (!with_cigar) { // score only + for (t = st_; t <= en_; ++t) { + __m256i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __dp_code_block1_pcl; + + z = _mm256_max_epi8(z, a); + z = _mm256_max_epi8(z, b); + z = _mm256_max_epi8(z, a2); + z = _mm256_max_epi8(z, b2); + z = _mm256_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; // save u[] and v[]; update a, b, a2 and b2 + if (t == en_) { + int ind = (en0 - t*32 + 16)/16;//doubt + // fprintf(stderr, "en0: %d, t: %d, ind: %d, msk: %d\n", en0, t, ind, msk); + //_mm256_storeu_si256(&u[t], (get_mask_store(msk_ar2_v[ind], &u[t], _mm256_sub_epi8(z, vt1)))); + _mm256_maskstore_epi32((int32_t *)&u[t], msk_ar2_v[ind], _mm256_sub_epi8(z, vt1)); + //_mm256_storeu_si256(&v[t], (get_mask_store(msk_ar2_v[ind], &v[t], _mm256_sub_epi8(z, ut)))); + _mm256_maskstore_epi32((int32_t *)&v[t], msk_ar2_v[ind], _mm256_sub_epi8(z, ut)); + tmp = _mm256_sub_epi8(z, q_); + a = _mm256_sub_epi8(a, tmp); + b = _mm256_sub_epi8(b, tmp); + tmp = _mm256_sub_epi8(z, q2_); + a2= _mm256_sub_epi8(a2, tmp); + b2= _mm256_sub_epi8(b2, tmp); + + } + else { + _mm256_storeu_si256(&u[t], _mm256_sub_epi8(z, vt1)); + _mm256_storeu_si256(&v[t], _mm256_sub_epi8(z, ut)); + tmp = _mm256_sub_epi8(z, q_); + a = _mm256_sub_epi8(a, tmp); + b = _mm256_sub_epi8(b, tmp); + tmp = _mm256_sub_epi8(z, q2_); + a2= _mm256_sub_epi8(a2, tmp); + b2= _mm256_sub_epi8(b2, tmp); + } + + if (t == en_) { + //__mmask32 msk; + int ind = (en0 - t*32 + 16)/16;//doubt + //assert(ind >= 0); + //msk = msk_ar2[ind]; + // fprintf(stderr, "en0: %d, t: %d, ind: %d, msk: %d\n", en0, t, ind, msk); + //_mm256_storeu_si256(&x[t], (get_mask_store(msk_ar2_v[ind], &x[t], _mm256_sub_epi8(_mm256_max_epi8(a, zero_), qe_)))); + _mm256_maskstore_epi32((int32_t *)&x[t], msk_ar2_v[ind], _mm256_sub_epi8(_mm256_max_epi8(a, zero_), qe_)); + //_mm256_storeu_si256(&y[t], (get_mask_store(msk_ar2_v[ind],&y[t], _mm256_sub_epi8(_mm256_max_epi8(b, zero_), qe_)))); + _mm256_maskstore_epi32((int32_t *)&y[t], msk_ar2_v[ind], _mm256_sub_epi8(_mm256_max_epi8(b, zero_), qe_)); + //_mm256_storeu_si256(&x2[t], (get_mask_store(msk_ar2_v[ind], &x2[t] , _mm256_sub_epi8(_mm256_max_epi8(a2, zero_), qe2_)))); + _mm256_maskstore_epi32((int32_t *)&x2[t], msk_ar2_v[ind], _mm256_sub_epi8(_mm256_max_epi8(a2, zero_), qe2_)); + //_mm256_storeu_si256(&y2[t], (get_mask_store(msk_ar2_v[ind], &y2[t], _mm256_sub_epi8(_mm256_max_epi8(b2, zero_), qe2_)))); + _mm256_maskstore_epi32((int32_t *)&y2[t], msk_ar2_v[ind], _mm256_sub_epi8(_mm256_max_epi8(b2, zero_), qe2_)); + } + else + { + _mm256_storeu_si256(&x[t], _mm256_sub_epi8(_mm256_max_epi8(a, zero_), qe_)); + _mm256_storeu_si256(&y[t], _mm256_sub_epi8(_mm256_max_epi8(b, zero_), qe_)); + _mm256_storeu_si256(&x2[t], _mm256_sub_epi8(_mm256_max_epi8(a2, zero_), qe2_)); + _mm256_storeu_si256(&y2[t], _mm256_sub_epi8(_mm256_max_epi8(b2, zero_), qe2_)); + } + // for (int l=0; l<64; l++) + // fprintf(stderr, "%d ", ((int8_t*)x)[l]); + + } + + } else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment + __m256i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + + + _mm_prefetch(pr + 4, _MM_HINT_T0); + for (t = st_; t < en_; ++t) { + __m256i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + // __mmask32 tmp_mask; + __m256i tmp_mask; + __dp_code_block1_pcl; + + d = _mm256_blendv_epi8(zero_, one_, _mm256_cmpgt_epi8(a, z)); // d = a > z? 1 : 0 + z = _mm256_max_epi8(z, a); + d = _mm256_blendv_epi8(d, two_,_mm256_cmpgt_epi8(b, z)); // d = b > z? 2 : d + z = _mm256_max_epi8(z, b); + d = _mm256_blendv_epi8(d, three_,_mm256_cmpgt_epi8(a2, z)); // d = a2 > z? 3 : d + z = _mm256_max_epi8(z, a2); + d = _mm256_blendv_epi8(d, four_,_mm256_cmpgt_epi8(b2, z)); // d = b2 > z? 4 : d + z = _mm256_max_epi8(z, b2); + z = _mm256_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + _mm256_storeu_si256(&u[t], _mm256_sub_epi8(z, vt1)); + _mm256_storeu_si256(&v[t], _mm256_sub_epi8(z, ut)); + tmp = _mm256_sub_epi8(z, q_); + a = _mm256_sub_epi8(a, tmp); + b = _mm256_sub_epi8(b, tmp); + tmp = _mm256_sub_epi8(z, q2_); + a2= _mm256_sub_epi8(a2, tmp); + b2= _mm256_sub_epi8(b2, tmp); + tmp_mask = _mm256_cmpgt_epi8(a, zero_); + _mm256_storeu_si256(&x[t], _mm256_sub_epi8(_mm256_blendv_epi8(zero_, a, tmp_mask), qe_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8( zero_, s1_, tmp_mask)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm256_cmpgt_epi8(b, zero_); + _mm256_storeu_si256(&y[t], _mm256_sub_epi8(_mm256_blendv_epi8( zero_, b, tmp_mask), qe_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8(zero_, s2_, tmp_mask)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm256_cmpgt_epi8(a2, zero_); + _mm256_storeu_si256(&x2[t], _mm256_sub_epi8(_mm256_blendv_epi8(zero_, a2,tmp_mask), qe2_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8( zero_, s3_,tmp_mask)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm256_cmpgt_epi8(b2, zero_); + _mm256_storeu_si256(&y2[t], _mm256_sub_epi8(_mm256_blendv_epi8(zero_, b2,tmp_mask), qe2_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8(zero_, s4_,tmp_mask)); // d = b > 0? 1<<6 : 0 + _mm256_storeu_si256(&pr[t], d); + _mm_prefetch(pr + 4, _MM_HINT_T0); + } + { + __m256i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + //__mmask32 tmp_mask; + __m256i tmp_mask; + __dp_code_block1_pcl; + + d = _mm256_blendv_epi8(zero_,one_, _mm256_cmpgt_epi8(a, z)); // d = a > z? 1 : 0 + z = _mm256_max_epi8(z, a); + d = _mm256_blendv_epi8(d, two_, _mm256_cmpgt_epi8(b, z)); // d = b > z? 2 : d + z = _mm256_max_epi8(z, b); + d = _mm256_blendv_epi8(d, three_, _mm256_cmpgt_epi8(a2, z)); // d = a2 > z? 3 : d + z = _mm256_max_epi8(z, a2); + d = _mm256_blendv_epi8(d, four_, _mm256_cmpgt_epi8(b2, z)); // d = b2 > z? 3 : d + z = _mm256_max_epi8(z, b2); + z = _mm256_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + { + //__mmask32 msk; + int ind = (en0 - t*32 + 16)/16;//doubt + // //assert(ind >= 0 && ind < 5); + //msk = msk_ar2[ind]; + + //_mm256_storeu_si256(&u[t], (get_mask_store(msk_ar2_v[ind], &u[t], _mm256_sub_epi8(z, vt1)))); + //_mm256_storeu_si256(&v[t], (get_mask_store(msk_ar2_v[ind],&v[t], _mm256_sub_epi8(z, ut)))); + _mm256_maskstore_epi32((int32_t *)&u[t], msk_ar2_v[ind], _mm256_sub_epi8(z, vt1)); + _mm256_maskstore_epi32((int32_t *)&v[t], msk_ar2_v[ind], _mm256_sub_epi8(z, ut)); + + + tmp = _mm256_sub_epi8(z, q_); + a = _mm256_sub_epi8(a, tmp); + b = _mm256_sub_epi8(b, tmp); + tmp = _mm256_sub_epi8(z, q2_); + a2= _mm256_sub_epi8(a2, tmp); + b2= _mm256_sub_epi8(b2, tmp); + } + + { + //__mmask32 msk; + int ind = (en0 - t*32 + 16)/16;//doubt + //msk = msk_ar2[ind]; + __m256i msk_v= msk_ar2_v[ind]; + off_end[r] -= (2-ind)*16;//doubt + + tmp_mask = _mm256_cmpgt_epi8(a, zero_); + //_mm256_storeu_si256(&x[t], (get_mask_store(msk_v, &x[t], _mm256_sub_epi8(_mm256_blendv_epi8(zero_, a, tmp_mask), qe_)))); + _mm256_maskstore_epi32((int32_t *)&x[t], msk_v, _mm256_sub_epi8(_mm256_blendv_epi8(zero_, a, tmp_mask), qe_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8(zero_, s1_, tmp_mask)); // d = a > 0? 1<<3 : 0 + + tmp_mask = _mm256_cmpgt_epi8(b, zero_); + //_mm256_storeu_si256(&y[t], (get_mask_store(msk_v,&y[t] , _mm256_sub_epi8(_mm256_blendv_epi8(zero_, b, tmp_mask), qe_)))); + _mm256_maskstore_epi32((int32_t *)&y[t], msk_v, _mm256_sub_epi8(_mm256_blendv_epi8(zero_, b, tmp_mask), qe_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8(zero_, s2_, tmp_mask)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm256_cmpgt_epi8(a2, zero_); + //_mm256_storeu_si256(&x2[t], (get_mask_store(msk_v,&x2[t], _mm256_sub_epi8(_mm256_blendv_epi8(zero_, a2, tmp_mask), qe2_)))); + _mm256_maskstore_epi32((int32_t *)&x2[t], msk_v, _mm256_sub_epi8(_mm256_blendv_epi8(zero_, a2, tmp_mask), qe2_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8(zero_, s3_, tmp_mask)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm256_cmpgt_epi8(b2, zero_); + //_mm256_storeu_si256(&y2[t], (get_mask_store(msk_v,&y2[t], _mm256_sub_epi8(_mm256_blendv_epi8(zero_, b2, tmp_mask), qe2_)))); + _mm256_maskstore_epi32((int32_t *)&y2[t], msk_v, _mm256_sub_epi8(_mm256_blendv_epi8(zero_, b2, tmp_mask), qe2_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8(zero_, s4_, tmp_mask)); // d = b > 0? 1<<6 : 0 + //_mm256_storeu_si256(&pr[t], (get_mask_store(msk_v, &pr[t], d))); + _mm256_maskstore_epi32((int32_t *)&pr[t], msk_v, d); + //_mm256_mask_storeu_epi8(&pr[t], msk, d); + } + } + + + } else { // gap right-alignment + __m256i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + // off[r] = stn, off_end[r] = enn; + // fprintf(stderr, "t: %d, st0: %d\n", st_, st0); + for (t = st_; t < en_; ++t) { + __m256i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + //__mmask32 tmp_mask; + __m256i tmp_mask; + + __dp_code_block1_pcl; + + d = _mm256_blendv_epi8(one_, zero_, _mm256_cmpgt_epi8(z, a)); + z = _mm256_max_epi8(z, a); + // d = _mm256_blendv_epi8(_mm256_cmpgt_epi8(z, b), _mm256_set1_epi8(2), d); + // d = z > b? d : 2 + d = _mm256_blendv_epi8(two_, d,_mm256_cmpgt_epi8(z, b) ); // d = z > b? d : 2 + z = _mm256_max_epi8(z, b); + // d = z > a2? d : 3 + d = _mm256_blendv_epi8(three_, d, _mm256_cmpgt_epi8(z, a2)); // d = z > a2? d : 3 + z = _mm256_max_epi8(z, a2); + // d = z > b2? d : 4 + d = _mm256_blendv_epi8(four_, d, _mm256_cmpgt_epi8(z, b2)); // d = z > b2? d : 4 + z = _mm256_max_epi8(z, b2); + z = _mm256_min_epi8(z, sc_mch_); + + // __dp_code_block2_pcl; + + //__mmask32 msk; + { + _mm256_storeu_si256(&u[t], _mm256_sub_epi8(z, vt1)); + _mm256_storeu_si256(&v[t], _mm256_sub_epi8(z, ut)); + tmp = _mm256_sub_epi8(z, q_); + a = _mm256_sub_epi8(a, tmp); + b = _mm256_sub_epi8(b, tmp); + tmp = _mm256_sub_epi8(z, q2_); + a2= _mm256_sub_epi8(a2, tmp); + b2= _mm256_sub_epi8(b2, tmp); + } + + { + tmp_mask = _mm256_cmpgt_epi8(zero_, a); + _mm256_storeu_si256(&x[t], _mm256_sub_epi8(_mm256_blendv_epi8(a, zero_,tmp_mask), qe_)); + // d = a > 0? 1<<3 : 0 + d = _mm256_or_si256(d, _mm256_blendv_epi8(s1_, zero_,tmp_mask)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm256_cmpgt_epi8(zero_, b); + _mm256_storeu_si256(&y[t], _mm256_sub_epi8(_mm256_blendv_epi8(b, zero_,tmp_mask), qe_)); + // d = b > 0? 1<<4 : 0 + d = _mm256_or_si256(d, _mm256_blendv_epi8(s2_, zero_,tmp_mask)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm256_cmpgt_epi8(zero_, a2); + _mm256_storeu_si256(&x2[t], _mm256_sub_epi8(_mm256_blendv_epi8(a2, zero_,tmp_mask), qe2_)); + // d = a > 0? 1<<5 : 0 + d = _mm256_or_si256(d, _mm256_blendv_epi8(s3_, zero_,tmp_mask)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm256_cmpgt_epi8(zero_, b2); + _mm256_storeu_si256(&y2[t], _mm256_sub_epi8(_mm256_blendv_epi8(b2, zero_,tmp_mask), qe2_)); + // d = b > 0? 1<<6 : 0 + d = _mm256_or_si256(d, _mm256_blendv_epi8(s4_, zero_,tmp_mask)); // d = b > 0? 1<<6 : 0 + _mm256_storeu_si256(&pr[t], d); + + } + } + //for (t = st_; t <= en_; ++t)// Last iteration unrolled + { + __m256i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + //__mmask32 tmp_mask; + __m256i tmp_mask; + + __dp_code_block1_pcl; + + d = _mm256_blendv_epi8(one_, zero_,_mm256_cmpgt_epi8(z, a) ); + z = _mm256_max_epi8(z, a); + // d = z > b? d : 2 + d = _mm256_blendv_epi8( two_, d, _mm256_cmpgt_epi8(z, b)); // d = z > b? d : 2 + z = _mm256_max_epi8(z, b); + // d = z > a2? d : 3 + d = _mm256_blendv_epi8( three_, d, _mm256_cmpgt_epi8(z, a2)); // d = z > a2? d : 3 + z = _mm256_max_epi8(z, a2); + // d = z > b2? d : 4 + d = _mm256_blendv_epi8(four_, d, _mm256_cmpgt_epi8(z, b2)); // d = z > b2? d : 4 + z = _mm256_max_epi8(z, b2); + z = _mm256_min_epi8(z, sc_mch_); + + // __dp_code_block2_pcl; + + //__mmask32 msk; + + // __mmask64 msk; + int ind = (en0 - t*32 + 16)/16;//doubt + //msk = msk_ar2[ind]; + off_end[r] -= (2-ind)*16;//doubt + + //_mm256_storeu_si256(&u[t], (get_mask_store(msk_ar2_v[ind], &u[t], _mm256_sub_epi8(z, vt1)))); + //_mm256_storeu_si256(&v[t], (get_mask_store(msk_ar2_v[ind], &v[t], _mm256_sub_epi8(z, ut)))); + _mm256_maskstore_epi32((int32_t *)&u[t], msk_ar2_v[ind], _mm256_sub_epi8(z, vt1)); + _mm256_maskstore_epi32((int32_t *)&v[t], msk_ar2_v[ind], _mm256_sub_epi8(z, ut)); + + + tmp = _mm256_sub_epi8(z, q_); + a = _mm256_sub_epi8(a, tmp); + b = _mm256_sub_epi8(b, tmp); + tmp = _mm256_sub_epi8(z, q2_); + a2= _mm256_sub_epi8(a2, tmp); + b2= _mm256_sub_epi8(b2, tmp); + + + + + tmp_mask = _mm256_cmpgt_epi8(zero_, a); + //_mm256_storeu_si256(&x[t], (get_mask_store(msk_ar2_v[ind],&x[t], _mm256_sub_epi8(_mm256_blendv_epi8(a, zero_, tmp_mask), qe_)))); + _mm256_maskstore_epi32((int32_t *)&x[t], msk_ar2_v[ind], _mm256_sub_epi8(_mm256_blendv_epi8(a, zero_,tmp_mask), qe_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8( s1_, zero_, tmp_mask)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm256_cmpgt_epi8(zero_, b); + //_mm256_storeu_si256(&y[t], (get_mask_store(msk_ar2_v[ind], &y[t], _mm256_sub_epi8(_mm256_blendv_epi8(b, zero_,tmp_mask), qe_)))); + _mm256_maskstore_epi32((int32_t *)&y[t], msk_ar2_v[ind], _mm256_sub_epi8(_mm256_blendv_epi8(b, zero_,tmp_mask), qe_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8(s2_, zero_,tmp_mask)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm256_cmpgt_epi8(zero_, a2); + //_mm256_storeu_si256(&x2[t], (get_mask_store(msk_ar2_v[ind], &x2[t], _mm256_sub_epi8(_mm256_blendv_epi8(a2, zero_,tmp_mask), qe2_)))); + _mm256_maskstore_epi32((int32_t *)&x2[t], msk_ar2_v[ind], _mm256_sub_epi8(_mm256_blendv_epi8(a2, zero_,tmp_mask), qe2_)); + d = _mm256_or_si256(d, _mm256_blendv_epi8(s3_, zero_,tmp_mask)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm256_cmpgt_epi8(zero_, b2); + //_mm256_storeu_si256(&y2[t], (get_mask_store(msk_ar2_v[ind], &y2[t], _mm256_sub_epi8(_mm256_blendv_epi8(b2, zero_,tmp_mask), qe2_)))); + _mm256_maskstore_epi32((int32_t *)&y2[t], msk_ar2_v[ind], _mm256_sub_epi8(_mm256_blendv_epi8(b2, zero_,tmp_mask), qe2_)); + + d = _mm256_or_si256(d, _mm256_blendv_epi8(s4_, zero_, tmp_mask)); // d = b > 0? 1<<6 : 0 + // _mm256_storeu_si256(&pr[t], d); +// _mm256_storeu_si256(&pr[t], (get_mask_store(msk_ar2_v[ind], &pr[t], d))); + _mm256_maskstore_epi32((int32_t *)&pr[t], msk_ar2_v[ind], d); + + + } + } + + if (!approx_max) { // find the exact max with a 32-bit score array + int32_t max_H, max_t; + // compute H[], max_H and max_t + if (r > 0) { + int32_t HH[8], tt[8], en1 = st0 + (en0 - st0) / 8 * 8, i; //doubt + __m256i max_H_, max_t_; + max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] : H[en0] + v8[en0]; // special casing the last element + + max_t = en0; + max_H_ = _mm256_set1_epi32(max_H); + max_t_ = _mm256_set1_epi32(max_t); +#if 1 + //for (t = st0; t < en1; t +=1) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t; + // H[t] += v8[t]; + //} + for (t = st0; t < en1; t += /*4*/8) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t; + __m256i H1, t_; + H1 = _mm256_loadu_si256((__m256i*)&H[t]); + __m128i t__ = _mm_load_si128((__m128i*) &v8[t]); + t_ = _mm256_cvtepi8_epi32(t__); + // t_ = _mm256_setr_epi32(v8[t], v8[t+1], v8[t+2], v8[t+3], v8[t+4], v8[t+5], v8[t+6], v8[t+7]); + H1 = _mm256_add_epi32(H1, t_); + _mm256_storeu_si256((__m256i*)&H[t], H1); + // making it 4 lanes to match accuracy + + __m256i shfH, shft, max1, max2; + __m256i tmp_mask_v; + t_ = _mm256_set1_epi32(t); + t_ = _mm256_add_epi32(t_, bt32_); + //shfH = _mm256_shuffle_i32x4(H1, H1, 0x1);//doubt + //shft = _mm256_shuffle_i32x4(t_, t_, 0x1);//doubt + shfH = _mm256_permute2x128_si256(H1, H1, 1); + shft = _mm256_permute2x128_si256(t_, t_, 1); + tmp_mask_v = _mm256_cmpgt_epi32(shfH, H1); + //max1 = _mm256_or_si256(_mm256_and_si256(tmp_mask_v, shfH) , _mm256_andnot_si256(tmp_mask_v, H1)); + max1 = _mm256_blendv_epi8(H1, shfH, tmp_mask_v); + //max2 = _mm256_or_si256(_mm256_and_si256(tmp_mask_v, shft) , _mm256_andnot_si256(tmp_mask_v, t_)); + max2 = _mm256_blendv_epi8(t_, shft, tmp_mask_v); + +/* //--shfH = _mm256_shuffle_i32x4(max1, max1, 0x2);//doubt + //--shft = _mm256_shuffle_i32x4(max2, max2, 0x2);//doubt + //--tmp_mask = _mm256_cmpgt_epi32_mask(shfH, max1); + //-max1 = _mm256_mask_blend_epi32(tmp_mask, max1, shfH); + //--max2 = _mm256_mask_blend_epi32(tmp_mask, max2, shft); + +*/ +// tmp_mask = _mm256_cmpgt_epi32_mask(max1, max_H_); +// max_H_ = _mm256_mask_blend_epi32(tmp_mask, max_H_, max1); +// max_t_ = _mm256_mask_blend_epi32(tmp_mask, max_t_, max2); + tmp_mask_v = _mm256_cmpgt_epi32(max1, max_H_); + //max_H_ = _mm256_or_si256(_mm256_and_si256(tmp_mask_v, max1) , _mm256_andnot_si256(tmp_mask_v, max_H_)); + //max_t_ = _mm256_or_si256(_mm256_and_si256(tmp_mask_v, max2) , _mm256_andnot_si256(tmp_mask_v, max_t_)); + max_H_ = _mm256_blendv_epi8(max_H_, max1, tmp_mask_v); + max_t_ = _mm256_blendv_epi8(max_t_, max2, tmp_mask_v); + + + } + _mm256_storeu_si256((__m256i*)HH, max_H_); + _mm256_storeu_si256((__m256i*)tt, max_t_); + + int rem = (en0 - t) / 4; + for (int l=0; l HH[j]) { + HH[j] = H[t]; + tt[j] = bt; + } + t++; + } + } + + for (i = 0; i < 4; ++i) + if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i; + for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE + H[t] += (int32_t)v8[t]; + if (H[t] > max_H) { + max_H = H[t], max_t = t; + } + } + #else + for (t = st0 ; t < en0; ++t) { // for the rest of values that haven't been computed with SSE + int32_t tmp = H[t] + v8[t]; + H[t] = tmp; + if (tmp > max_H) { + max_H = tmp, max_t = t; + } + } + #endif + } else H[0] = v8[0] - qe, max_H = H[0], max_t = 0; // special casing r==0 + + // update ez + if (en0 == tlen - 1 && H[en0] > ez->mte) { + ez->mte = H[en0], ez->mte_q = r - en_new; + } + if (r - st0 == qlen - 1 && H[st0] > ez->mqe) { + ez->mqe = H[st0], ez->mqe_t = st0; + } + + if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e2)) { + break; + } + if (r == qlen + tlen - 2 && en0 == tlen - 1) { + ez->score = H[tlen - 1]; + } + } else { // find approximate max; Z-drop might be inaccurate, too. + if (r > 0) { + if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) { + int32_t d0 = v8[last_H0_t]; + int32_t d1 = u8[last_H0_t + 1]; + if (d0 > d1) H0 += d0; + else H0 += d1, ++last_H0_t; + } else if (last_H0_t >= st0 && last_H0_t <= en0) { + H0 += v8[last_H0_t]; + } else { + ++last_H0_t, H0 += u8[last_H0_t]; + } + } else H0 = v8[0] - qe, last_H0_t = 0; + if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e2)) { + break; + } + if (r == qlen + tlen - 2 && en0 == tlen - 1) { + ez->score = H0; + } + } + // last_st = st, last_en = en; + last_st = st_new, last_en = en_new; + + } +#ifdef MANUAL_PROFILING + alignment_time += (__rdtsc() - align_start); +#endif + + kfree(km, mem); + if (!approx_max) kfree(km, H); + if (with_cigar) { // backtrack + int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR); + if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*32, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) { + ez->reach_end = 1; + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*32, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (ez->max_t >= 0 && ez->max_q >= 0) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*32, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } + kfree(km, mem2); kfree(km, off); + } +// kfree_all(); +#undef __dp_code_block1_pcl +#undef __dp_code_block2_pcl +} + +#endif + diff --git a/ksw2_extd2_avx.h b/ksw2_extd2_avx.h new file mode 100644 index 0000000..69db90b --- /dev/null +++ b/ksw2_extd2_avx.h @@ -0,0 +1,42 @@ +/* The MIT License + +Copyright (c) 2018- Dana-Farber Cancer Institute + 2017-2018 Broad Institute, Inc. + +Permission is hereby granted, free of charge, to any person obtaining +a copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be +included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS +BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN +ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +Modified Copyright (C) 2021 Intel Corporation + Contacts: Saurabh Kalikar ; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ +#include +#include +#include +#include "ksw2.h" +#include +#include +#include +#include +void ksw_extd2_avx512(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez); + +void ksw_extd2_avx2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez); diff --git a/lchain.c b/lchain.c index a1615ac..febd636 100644 --- a/lchain.c +++ b/lchain.c @@ -5,6 +5,11 @@ #include "mmpriv.h" #include "kalloc.h" #include "krmq.h" +//#include "simd_chain.h" +#include "parallel_chaining_32_bit.h" + +extern uint64_t dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4; + uint64_t *mg_chain_backtrack(void *km, int64_t n, const int32_t *f, const int64_t *p, int32_t *v, int32_t *t, int32_t min_cnt, int32_t min_sc, int32_t *n_u_, int32_t *n_v_) { @@ -88,15 +93,65 @@ static mm128_t *compact_a(void *km, int32_t n_u, uint64_t *u, int32_t n_v, int32 static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t max_dist_x, int32_t max_dist_y, int32_t bw, float chn_pen_gap, float chn_pen_skip, int is_cdna, int n_seg) { - int32_t dq = (int32_t)ai->y - (int32_t)aj->y, dr, dd, dg, q_span, sc; - int32_t sidi = (ai->y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT; - int32_t sidj = (aj->y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT; - if (dq <= 0 || dq > max_dist_x) return INT32_MIN; - dr = (int32_t)(ai->x - aj->x); - if (sidi == sidj && (dr == 0 || dq > max_dist_y)) return INT32_MIN; + + uint64_t ai_x, ai_y, aj_x, aj_y; + ai_x = ai->x; ai_y = ai->y; aj_x = aj->x; aj_y = aj->y; + +#ifdef CHAIN_DEBUG + int32_t sc_vect = obj.comput_sc_vectorized_avx2_caller(ai_x, ai_y, aj_x, aj_y, aj->y>>32&0xff); +#endif + + //if (sc_vect == 0) return INT32_MIN; + //else + //return sc_vect; + + //fprintf(stderr, "%lld %lld %lld %lld \n", ai_x, ai_y, aj_x, aj_y); + //fprintf(stderr, "%lld %lld %lld %f %f %d %d\n", max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg); + int32_t dq = (int32_t)ai_y - (int32_t)aj_y, dr, dd, dg, q_span, sc; + int32_t sidi = (ai_y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT; + int32_t sidj = (aj_y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT; + if (dq <= 0 || dq > max_dist_x) { + +#ifdef CHAIN_DEBUG + if(INT32_MIN != sc_vect){ + //fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect); + fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect); + } +#endif + return INT32_MIN; + } + dr = (int32_t)(ai_x - aj_x); + if (sidi == sidj && (dr == 0 || dq > max_dist_y)) { + +#ifdef CHAIN_DEBUG + if(INT32_MIN != sc_vect){ + //fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect); + fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect); + } +#endif + return INT32_MIN; + } dd = dr > dq? dr - dq : dq - dr; - if (sidi == sidj && dd > bw) return INT32_MIN; - if (n_seg > 1 && !is_cdna && sidi == sidj && dr > max_dist_y) return INT32_MIN; + if (sidi == sidj && dd > bw) { + +#ifdef CHAIN_DEBUG + if(INT32_MIN != sc_vect){ + //fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect); + fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect); + } +#endif + return INT32_MIN; + } + if (n_seg > 1 && !is_cdna && sidi == sidj && dr > max_dist_y) { + +#ifdef CHAIN_DEBUG + if(INT32_MIN != sc_vect){ + //fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect); + fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect); + } +#endif + return INT32_MIN; + } dg = dr < dq? dr : dq; q_span = aj->y>>32&0xff; sc = q_span < dg? q_span : dg; @@ -110,6 +165,13 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma else sc -= (int)(lin_pen + .5f * log_pen); } else sc -= (int)(lin_pen + .5f * log_pen); } +#ifdef CHAIN_DEBUG + + if(sc != sc_vect ){ + //fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect); + fprintf(stderr, "outer: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect); + } +#endif return sc; } @@ -124,10 +186,18 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip, int is_cdna, int n_seg, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km) { // TODO: make sure this works when n has more than 32 bits - int32_t *f, *t, *v, n_u, n_v, mmax_f = 0; + ///fprintf(stderr, "chaining called\n"); + + + +#ifdef MANUAL_PROFILING + uint64_t align_start = __rdtsc(); +#endif + + int32_t *f, *t, *v, *v_1, *p_1, n_u, n_v, mmax_f = 0; int64_t *p, i, j, max_ii, st = 0, n_iter = 0; uint64_t *u; - + uint32_t* f_1; if (_u) *_u = 0, *n_u_ = 0; if (n == 0 || a == 0) { kfree(km, a); @@ -136,20 +206,67 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int if (max_dist_x < bw) max_dist_x = bw; if (max_dist_y < bw && !is_cdna) max_dist_y = bw; KMALLOC(km, p, n); + KMALLOC(km, p_1, n); KMALLOC(km, f, n); + KMALLOC(km, f_1, n); KMALLOC(km, v, n); + KMALLOC(km, v_1, n); KCALLOC(km, t, n); +#ifdef PARALLEL_CHAINING + +// Parallel chaining data-structures + anchor_t* anchors = (anchor_t*)malloc(n* sizeof(anchor_t)); + for (i = 0; i < n; ++i) { + uint64_t ri = a[i].x; + int32_t qi = (int32_t)a[i].y, q_span = a[i].y>>32&0xff; // NB: only 8 bits of span is used!!! + anchors[i].r = ri; + anchors[i].q = qi; + anchors[i].l = q_span; + } + num_bits_t *anchor_r, *anchor_q, *anchor_l; + create_SoA_Anchors_32_bit(anchors, n, anchor_r, anchor_q, anchor_l); + //dp_chain obj(max_dist_x, max_dist_y, bw, max_skip, max_iter, 0, is_cdna, n_seg); + dp_chain obj(max_dist_x, max_dist_y, bw, max_skip, max_iter, min_cnt, min_sc, chn_pen_gap, chn_pen_skip, is_cdna, n_seg); + + obj.mm_dp_vectorized(n, &anchors[0], anchor_r, anchor_q, anchor_l, f_1, p_1, v_1, max_dist_x, max_dist_y, NULL, NULL); + + // -16 is due to extra padding at the start of arrays + anchor_r -= 16; anchor_q -= 16; anchor_l -= 16; + free(anchor_r); + free(anchor_q); + free(anchor_l); + free(anchors); + for(int i = 0; i < n; i++){ + // if(f[i] != f_1[i] || p[i] != p_1[i] || v[i] !=v_1[i]) + // { +// fprintf(stderr, "i:%d %d %d %d %d %d %d\n",i, f[i], f_1[i], p[i], p_1[i], v[i], v_1[i] ); + // } +#if 1 + f[i] = f_1[i]; + p[i] = p_1[i]; + v[i] = v_1[i]; +#endif + } + +// +#else + // fill the score and backtrack arrays for (i = 0, max_ii = -1; i < n; ++i) { int64_t max_j = -1, end_j; int32_t max_f = a[i].y>>32&0xff, n_skip = 0; while (st < i && (a[i].x>>32 != a[st].x>>32 || a[i].x > a[st].x + max_dist_x)) ++st; if (i - st > max_iter) st = i - max_iter; + int my_cnt = 0; for (j = i - 1; j >= st; --j) { int32_t sc; sc = comput_sc(&a[i], &a[j], max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg); ++n_iter; + // if(i == 177){ + //fprintf(stderr, "args: %d %d %d %d %d\n", a[i].x, a[i].y, a[j].x, a[j].y, a[j].y>>32&0xff); + //fprintf(stderr, "j_th %d score: %d\n", ++my_cnt, sc); + // } if (sc == INT32_MIN) continue; sc += f[j]; if (sc > max_f) { @@ -162,32 +279,72 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int if (p[j] >= 0) t[p[j]] = i; } end_j = j; + int debug_iter = 2057329; + //if (i == debug_iter) fprintf(stderr, "mm2 -- endj: %d max_ii: %d max_f: %d \n", end_j, max_ii, max_f); + +#if 1 if (max_ii < 0 || a[i].x - a[max_ii].x > (int64_t)max_dist_x) { int32_t max = INT32_MIN; max_ii = -1; - for (j = i - 1; j >= st; --j) - if (max < f[j]) max = f[j], max_ii = j; + for (j = i - 1; j >= st; --j) { + if (max < (int32_t)f[j]) max = f[j], max_ii = j; + } } +#endif +#if 1 if (max_ii >= 0 && max_ii < end_j) { int32_t tmp; tmp = comput_sc(&a[i], &a[max_ii], max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg); - if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) + // if (i == debug_iter) fprintf(stderr, "mm2: endj: %d max_ii: %d max_f: %d tmp_score: %d \n", end_j, max_ii, max_f, tmp); + + + if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) { + // if (i == debug_iter) fprintf(stderr, "mm2: endj: %d max_ii: %d max_f: %d tmp_score: %d \n", end_j, max_ii, max_f, tmp); max_f = tmp + f[max_ii], max_j = max_ii; + // if (i == debug_iter) fprintf(stderr, "mm2: endj: %d max_ii: %d max_f: %d tmp_score: %d sum : %d \n", end_j, max_ii, max_f, tmp, tmp + f[max_ii]); + + } } +#endif f[i] = max_f, p[i] = max_j; v[i] = max_j >= 0 && v[max_j] > max_f? v[max_j] : max_f; // v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak + +#if 1 if (max_ii < 0 || (a[i].x - a[max_ii].x <= (int64_t)max_dist_x && f[max_ii] < f[i])) max_ii = i; if (mmax_f < max_f) mmax_f = max_f; +#endif } +#endif + +#ifdef CHAIN_DEBUG + + for(int i = 0; i < n; i++){ + if(f[i] != f_1[i] || p[i] != p_1[i] || v[i] !=v_1[i]) + { + fprintf(stderr, "i:%d %d %d %d %d %d %d\n",i, f[i], f_1[i], p[i], p_1[i], v[i], v_1[i] ); + } +#if 0 + f[i] = f_1[i]; + p[i] = p_1[i]; + v[i] = v_1[i]; +#endif + } + +#endif u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, &n_u, &n_v); *n_u_ = n_u, *_u = u; // NB: note that u[] may not be sorted by score here - kfree(km, p); kfree(km, f); kfree(km, t); + kfree(km, p); kfree(km, p_1); kfree(km, f); kfree(km, f_1); kfree(km, t); kfree(km, v_1); if (n_u == 0) { kfree(km, a); kfree(km, v); return 0; } + + +#ifdef MANUAL_PROFILING + dp_time += __rdtsc() - align_start; +#endif return compact_a(km, n_u, u, n_v, v, a); } @@ -225,6 +382,11 @@ static inline int32_t comput_sc_simple(const mm128_t *ai, const mm128_t *aj, flo mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_skip, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km) { +#ifdef MANUAL_PROFILING + uint64_t start = __rdtsc(); +#endif + uint64_t tim; + //fprintf(stderr, "rmq call \n"); int32_t *f,*t, *v, n_u, n_v, mmax_f = 0, max_rmq_size = 0; int64_t *p, i, i0, st = 0, st_inner = 0, n_iter = 0; uint64_t *u; @@ -252,6 +414,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski int32_t q_span = a[i].y>>32&0xff, max_f = q_span; lc_elem_t s, *q, *r, lo, hi; // add in-range anchors +#ifdef MANUAL_PROFILING_RMQ + tim = __rdtsc(); +#endif if (i0 < i && a[i0].x != a[i].x) { int64_t j; for (j = i0; j < i; ++j) { @@ -266,7 +431,13 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski } i0 = i; } +#ifdef MANUAL_PROFILING_RMQ + rmq_t1 += __rdtsc() - tim; +#endif // get rid of active chains out of range +#ifdef MANUAL_PROFILING_RMQ + tim = __rdtsc(); +#endif while (st < i && (a[i].x>>32 != a[st].x>>32 || a[i].x > a[st].x + max_dist || krmq_size(head, root) > cap_rmq_size)) { s.y = (int32_t)a[st].y, s.i = st; if ((q = krmq_find(lc_elem, root, &s, 0)) != 0) { @@ -275,6 +446,12 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski } ++st; } +#ifdef MANUAL_PROFILING_RMQ + rmq_t2 += __rdtsc() - tim; +#endif +#ifdef MANUAL_PROFILING_RMQ + tim = __rdtsc(); +#endif if (max_dist_inner > 0) { // similar to the block above, but applied to the inner tree while (st_inner < i && (a[i].x>>32 != a[st_inner].x>>32 || a[i].x > a[st_inner].x + max_dist_inner || krmq_size(head, root_inner) > cap_rmq_size)) { s.y = (int32_t)a[st_inner].y, s.i = st_inner; @@ -285,6 +462,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski ++st_inner; } } +#ifdef MANUAL_PROFILING_RMQ + rmq_t3 += __rdtsc() - tim; +#endif // RMQ lo.i = INT32_MAX, lo.y = (int32_t)a[i].y - max_dist; hi.i = 0, hi.y = (int32_t)a[i].y; @@ -304,6 +484,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski krmq_itr_t(lc_elem) itr; krmq_itr_find(lc_elem, root_inner, lo, &itr); while ((q = krmq_at(&itr)) != 0) { +#ifdef MANUAL_PROFILING_RMQ + tim = __rdtsc(); +#endif if (q->y < (int32_t)a[i].y - max_dist_inner) break; ++n_rmq_iter; j = q->i; @@ -319,11 +502,15 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski if (p[j] >= 0) t[p[j]] = i; } if (!krmq_itr_prev(lc_elem, &itr)) break; +#ifdef MANUAL_PROFILING_RMQ + rmq_t4 += __rdtsc() - tim; +#endif } n_iter += n_rmq_iter; } } } + // set max assert(max_j < 0 || (a[max_j].x < a[i].x && (int32_t)a[max_j].y < (int32_t)a[i].y)); f[i] = max_f, p[i] = max_j; @@ -340,5 +527,8 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski kfree(km, a); kfree(km, v); return 0; } +#ifdef MANUAL_PROFILING + rmq_time += __rdtsc() - start; +#endif return compact_a(km, n_u, u, n_v, v, a); } diff --git a/lisa_hash.h b/lisa_hash.h new file mode 100644 index 0000000..5d1b89a --- /dev/null +++ b/lisa_hash.h @@ -0,0 +1,307 @@ +/************************************************************************************* +MIT License + +Copyright (c) 2020 Intel Labs + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + +Authors: Saurabh Kalikar ; Sanchit Misra ; +*****************************************************************************************/ + +#include +#include +#include +#include +#include +#include "RMI.h" +using namespace std; + +template +class lisa_hash{ + public: + + rmi_val_t* p; + rmi_val_t* p_bin; + + private: + + rmi_key_t* keys; + + uint64_t* values_enc; + uint64_t* values_enc_bin; + + uint64_t keys_size; + uint64_t p_size; + + RMI *rmi; + + void mem_alloc(uint64_t hash_size, uint64_t p_size){ + + // Removed for memory optimization + #if 0 + keys = (rmi_key_t*)malloc((1+hash_size)*sizeof(rmi_key_t)); + keys = keys + 1; // keys[-1] stores total number of keys + #endif + values_enc = (uint64_t*)malloc(hash_size*sizeof(uint64_t)); + + p = (rmi_val_t*)malloc(p_size*sizeof(rmi_val_t)); + + //values_enc_bin = (uint64_t*)malloc(hash_size*sizeof(uint64_t)); + + //p_bin = (rmi_val_t*)malloc(p_size*sizeof(rmi_val_t)); + + keys_size = hash_size; + // keys[-1] = keys_size; + this->p_size = p_size; + + fprintf(stderr, "Memory allocated %lld \n", p_size); + + } + //This function is no longer used. + void store_sorted_keys(string binFileName){ + + + ofstream wf(binFileName, ios::out | ios::binary); + wf.write((char*)&keys[-1], (1+keys_size)*sizeof(rmi_key_t)); + + wf.close(); + + rmi_key_t *temp_keys = (rmi_key_t*)malloc((1 + keys_size)*sizeof(rmi_key_t)); + ifstream rf(binFileName, ios::out | ios::binary); + rf.read((char*)&temp_keys[0], (1 + keys_size)*sizeof(rmi_key_t)); + + if(temp_keys[0] == keys[-1]){ + for(int i = 0; i < keys_size; i++){ + if(temp_keys[i+1] != keys[i]){ + cout<<"Error: "<> 32; + p_size = (uint32_t) values_enc[i]; + } + + + + void load_bin(string inputFile){ + fprintf(stderr, "Loading from bin\n"); + + string f1_name = (string) inputFile + "_pos_bin"; + string f2_name = (string) inputFile + "_val_bin"; + ifstream instream_f1(f1_name, ifstream::binary); + ifstream instream_f2(f2_name, ifstream::binary); + instream_f1.seekg(0); + instream_f1.read((char*)&values_enc[0], keys_size*sizeof(uint64_t)); + instream_f1.close(); + + instream_f2.seekg(0); + instream_f2.read((char*)&p[0], p_size*sizeof(uint64_t)); + instream_f2.close(); + /* + for(uint64_t i = 0; i < p_size; i++){ + if(p[i] != p_bin[i]){ + fprintf(stderr, "Error!! %lld %lld\n",p[i], p_bin[i] ); + } + } + for(uint64_t i = 0; i < keys_size; i++){ + if(values_enc[i] != values_enc_bin[i]){ + fprintf(stderr, "Error!! %lld %lld\n", values_enc[i], values_enc_bin[i]); + } + } + */ + } + + void load(string inputFile){ + + ifstream f(inputFile); + + rmi_key_t key; + rmi_val_t val; + uint64_t n; + int64_t offset = 0; + int64_t i = 0; + while(f>>key){ + f>>n; + //values_enc[++key_size] = val_size << 32 | n; + for(int j = 0; j < n; j++){ + f>>val; + p[offset + j] = val; + } + set_key_value(i, key, offset, n); + offset +=n; + i++; + } + } + + + + public: + lisa_hash(string inputFile, char* rmi_prefix, long leaf = 0){ + uint64_t start_time, load_time, rmi_building_time, rmi_object; + start_time = __rdtsc(); + + uint64_t val_count = 0; + + ifstream f_size(inputFile+"_size"); + f_size>>keys_size; + f_size>>p_size; + mem_alloc(keys_size, p_size); + + fprintf(stderr, "Num_keys: %lld, num_values = %lld", keys_size, p_size); + //load(inputFile); + load_bin(inputFile); + + string prefix = inputFile + "_keys"; +#ifdef UINT64 + string keys_bin_file_name = prefix + ".uint64"; +#else + string keys_bin_file_name = prefix + ".f64"; +#endif + + ifstream rf(keys_bin_file_name, ios::out | ios::binary); + if(!rf.good()){ + cout<<"Error: Binary file with keys not found!!\n"; + //call store keys + //store_sorted_keys(keys_bin_file_name); + exit(0); + } + + + string keys_rmi_file = prefix + ".rmi_PARAMETERS"; + ifstream rmi_f(keys_rmi_file, ios::out | ios::binary); + if(!rmi_f.good() || leaf != 0){ + cout<<"rmi file not found: "<< keys_rmi_file <(&prefix[0]); + + + } + + + rmi_val_t* get_hash_value(rmi_key_t key, int*n){ + uint64_t index = rmi->lookup(key); + if(index == -1){ + cout<<"Key not found\n"; + return NULL; + } + uint64_t pos, p_size; + + get_val(index, pos, p_size); + *n = p_size; + return &p[pos]; + } + + rmi_val_t* get_hash_values_batched(rmi_key_t *keys, uint64_t num_keys, int* &num_values){ + + int64_t *pos = (int64_t*) malloc(num_keys*sizeof(int64_t)); + + rmi->lookup_batched(keys, num_keys, &pos[0]); + + + num_values = (int*) malloc(num_keys*sizeof(int)); + rmi_val_t **p_ptrs = (rmi_val_t**) malloc(num_keys*sizeof(rmi_val_t*)); + uint64_t total_num_values = 0; + for(int i = 0; i < num_keys; i++){ + if(pos[i] == -1){ + num_values[i] = 0; + continue; + } + num_values[i] = (uint32_t) values_enc[pos[i]]; + total_num_values+= num_values[i]; + p_ptrs[i] = &p[values_enc[pos[i]] >> 32]; + } + free(pos); + rmi_val_t* ret_values = (rmi_val_t*) malloc(total_num_values*sizeof(rmi_val_t)); + + uint64_t cnt= 0; + for(int i = 0; i < num_keys; i++){ + + rmi_val_t *value_start_ptr = p_ptrs[i]; + uint64_t numhit = num_values[i]; + + + for(int j = 0; j < numhit; j++){ + ret_values[cnt++] = value_start_ptr[j]; + } + } + free(p_ptrs); + return ret_values; + } + + + void mm_idx_get_batched(uint64_t* &minimizers, uint64_t num_minimizers, int64_t* &pos, uint64_t** &p_ptrs, int* &num_hits){ + + rmi->lookup_batched(minimizers, num_minimizers, &pos[0]); + + for(int i = 0; i < num_minimizers; i++){ + int64_t p_i = pos[i]; + + if(p_i < 0 || p_i > keys_size ) + num_hits[i] = 0; + else + num_hits[i] = (uint32_t) values_enc[p_i]; + p_ptrs[i] = p + (values_enc[p_i] >> 32); + } + } + ~lisa_hash(){ + delete rmi; + free(values_enc); + free(p); + } + +}; diff --git a/main.c b/main.c index 559971e..0554d84 100644 --- a/main.c +++ b/main.c @@ -1,5 +1,3 @@ -/* Fast contrib: v.22 */ - #include #include #include @@ -8,6 +6,73 @@ #include "minimap.h" #include "mmpriv.h" #include "ketopt.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "bseq.h" +#include "minimap.h" +#include "mmpriv.h" +#include "ketopt.h" + +//#include "profile.h" +#include +#include +#include + +using namespace std; +uint64_t avg; +uint64_t minimizer_lookup_time, alignment_time, dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4; + +#ifdef LISA_HASH +#include "lisa_hash.h" +lisa_hash *lh; +#endif + +// New memory allocation approach for alignment optimizations +// +void *km1; +uint64_t km_size = 500000000; // 500 MB +int km_top; +/* +void *kcalloc_(void* km, int count, int size) +{ + assert(count*size < km_size); + km_top += count*size + 1024; + memset(km, 0, count * size); + + // printf("km_top: %d\n", km_top); + return km; +} + +void *kmalloc_(void* km, int count) { + if(km_top + count >= km_size) + printf("count: %d\n", count); + assert(km_top + count < km_size); + void *mem = (void*) ((int8_t*) km + km_top); + km_top += count + 1024; + // printf("km_top: %d\n", km_top); + return mem; +} + +void kfree_all() { km_top = 0;} + */ + +// Memory for alignment end + + +#ifndef __rdtsc +#ifdef _rdtsc +#define __rdtsc _rdtsc +#else +#define __rdtsc __builtin_ia32_rdtsc +#endif +#endif #define MM_VERSION "2.22-r1101" @@ -119,6 +184,9 @@ static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const int main(int argc, char *argv[]) { + // Memory allocation for alignment optimizations + //km1 = calloc(km_size, 1); // 10 MB init contg. alloc + const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:"; ketopt_t o = KETOPT_INIT; mm_mapopt_t opt; @@ -133,9 +201,11 @@ int main(int argc, char *argv[]) liftrlimit(); mm_realtime0 = realtime(); mm_set_opt(0, &ipt, &opt); + string preset_arg = ""; while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { // test command line options and apply option -x/preset first if (c == 'x') { + preset_arg += (string) o.arg; if (mm_set_opt(o.arg, &ipt, &opt) < 0) { fprintf(stderr, "[ERROR] unknown preset '%s'\n", o.arg); return 1; @@ -368,6 +438,7 @@ int main(int argc, char *argv[]) fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n"); return 1; } + preset_arg = (string)argv[o.ind] + "_" + preset_arg + "_minimizers_key_value_sorted"; idx_rdr = mm_idx_reader_open(argv[o.ind], &ipt, fnw); if (idx_rdr == 0) { fprintf(stderr, "[ERROR] failed to open file '%s': %s\n", argv[o.ind], strerror(errno)); @@ -410,6 +481,9 @@ int main(int argc, char *argv[]) __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq); if (argc != o.ind + 1) mm_mapopt_update(&opt, mi); if (mm_verbose >= 3) mm_idx_stat(mi); +#ifdef LISA_INDEX + mm_idx_dump_hash(preset_arg.c_str(), mi); +#endif if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1); if (alt_list) mm_idx_alt_read(mi, alt_list); if (argc - (o.ind + 1) == 0) { @@ -417,6 +491,16 @@ int main(int argc, char *argv[]) continue; // no query files } ret = 0; +#ifdef LISA_HASH + fprintf(stderr, "Using LISA_HASH..\n"); + mm_idx_destroy_mm_hash(mi); + char* prefix; + lh = new lisa_hash(preset_arg, prefix); + fprintf(stderr, "Loading done.\n"); +// total_time = __rdtsc(); +// fprintf(stderr, "\nIndexing Real time: %.3f sec;\n", realtime() - mapping_time); +#endif + mm_realtime0 = realtime(); if (!(opt.flag & MM_F_FRAG_MODE)) { for (i = o.ind + 1; i < argc; ++i) { ret = mm_map_file(mi, argv[i], &opt, n_threads); @@ -425,12 +509,17 @@ int main(int argc, char *argv[]) } else { ret = mm_map_file_frag(mi, argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_threads); } - mm_idx_destroy(mi); + //mm_idx_destroy(mi); if (ret < 0) { fprintf(stderr, "ERROR: failed to map the query file\n"); exit(EXIT_FAILURE); } } +#ifdef LISA_HASH + mm_idx_destroy_seq(mi); +#else + mm_idx_destroy(mi); +#endif n_parts = idx_rdr->n_parts; mm_idx_reader_close(idx_rdr); @@ -449,5 +538,10 @@ int main(int argc, char *argv[]) fprintf(stderr, " %s", argv[i]); fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0); } + + fprintf(stderr, "minimizer-lookup: %lld dp: %lld rmq: %lld rmq_t1: %lld rmq_t2: %lld rmq_t3: %lld rmq_t4: %lld alignment: %lld %lld\n", minimizer_lookup_time, dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4, alignment_time, avg); +#ifdef LISA_HASH + delete lh; +#endif return 0; } diff --git a/map.c b/map.c index 9c4ad64..238bd6b 100644 --- a/map.c +++ b/map.c @@ -10,6 +10,8 @@ #include "bseq.h" #include "khash.h" + +extern uint64_t rmq_time; struct mm_tbuf_s { void *km; int rep_len, frag_gap; @@ -276,7 +278,12 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** if (opt->flag & MM_F_RMQ) { a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km); +// a = mg_lchain_dp(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score, +// opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); + + } else { + //fprintf(stderr, "dp call - n_a = %lld\n", n_a); a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); } @@ -284,12 +291,23 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** if (opt->bw_long > opt->bw && (opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1) { // re-chain/long-join for long sequences int32_t st = (int32_t)a[0].y, en = (int32_t)a[(int32_t)u[0] - 1].y; if (qlen_sum - (en - st) > opt->rmq_rescue_size || en - st > qlen_sum * opt->rmq_rescue_ratio) { +#ifdef MANUAL_PROFILING +// uint64_t tim = __rdtsc(); +#endif + // fprintf(stderr, "pre: rmq rechain call - n_a = %lld n_regs = %lld\n",n_a, n_regs0); int32_t i; + int64_t prev_n_a = n_a; for (i = 0, n_a = 0; i < n_regs0; ++i) n_a += (int32_t)u[i]; kfree(b->km, u); radix_sort_128x(a, a + n_a); +// fprintf(stderr, "post: rmq rechain call - prev_n_a = %lld n_a = %lld n_regs = %lld\n",prev_n_a, n_a, n_regs0); +// a = mg_lchain_dp(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score, +// opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km); +#ifdef MANUAL_PROFILING +// rmq_time += __rdtsc() - tim; +#endif } } else if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ)) { // re-chain, mostly for short reads int rechain = 0; @@ -563,6 +581,7 @@ static void *worker_pipeline(void *shared, int step, void *in) if ((p->opt->flag & MM_F_OUT_CS) && !(mm_dbg_flag & MM_DBG_NO_KALLOC)) km = km_init(); for (k = 0; k < s->n_frag; ++k) { int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k]; +#ifndef DISABLE_OUTPUT for (i = seg_st; i < seg_en; ++i) { mm_bseq1_t *t = &s->seq[i]; if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files @@ -597,6 +616,7 @@ static void *worker_pipeline(void *shared, int step, void *in) mm_err_puts(p->str.s); } } +#endif for (i = seg_st; i < seg_en; ++i) { for (j = 0; j < s->n_reg[i]; ++j) free(s->reg[i][j].p); free(s->reg[i]); diff --git a/minimap.h b/minimap.h index a7a886b..a7f5d62 100644 --- a/minimap.h +++ b/minimap.h @@ -285,6 +285,13 @@ mm_idx_t *mm_idx_load(FILE *fp); */ void mm_idx_dump(FILE *fp, const mm_idx_t *mi); +/** + * Store hash table from minimap2 index into a file + * @param f_name File name for output file + * @param mi minimap2 index + */ +void mm_idx_dump_hash(const char* f_name, const mm_idx_t *mi); + /** * Create an index from strings in memory * @@ -313,6 +320,19 @@ void mm_idx_stat(const mm_idx_t *idx); * @param r minimap2 index */ void mm_idx_destroy(mm_idx_t *mi); +/** + * Destroy/deallocate an hash table index + * + * @param r minimap2 index + */ +void mm_idx_destroy_mm_hash(mm_idx_t *mi); + +/** + * Destroy/deallocate target sequences + * + * @param r minimap2 index + */ +void mm_idx_destroy_seq(mm_idx_t *mi); /** * Initialize a thread-local buffer for mapping diff --git a/seed.c b/seed.c index baf1b20..a595ef0 100644 --- a/seed.c +++ b/seed.c @@ -1,9 +1,40 @@ #include "mmpriv.h" #include "kalloc.h" #include "ksort.h" +#include +#include +#include + +#ifdef LISA_HASH +#include "lisa_hash.h" +extern lisa_hash *lh; +#endif +extern uint64_t minimizer_lookup_time; mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, int32_t *n_m_) { +#ifdef MANUAL_PROFILING + uint64_t lookup_start = __rdtsc(); +#endif + +#ifdef LISA_HASH +//----------------------------------- + uint64_t** cr_batch = (uint64_t**) malloc((mv->n)*sizeof(uint64_t*)); + int* t_batch = (int*)malloc((mv->n)*sizeof(int)); + uint64_t* minimizers = (uint64_t*) malloc((mv->n)*sizeof(uint64_t)); + int64_t* lisa_pos = (int64_t*) malloc((max(32, (int)mv->n))* sizeof(int64_t)); + + for (size_t i = 0; i < mv->n; i++) { + mm128_t *p = &mv->a[i]; + minimizers[i] = p->x>>8; + } + + lh->mm_idx_get_batched(minimizers, mv->n, lisa_pos, cr_batch, t_batch); +//----------------------------------- + +#endif + + mm_seed_t *m; size_t i; int32_t k; @@ -14,7 +45,12 @@ mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, mm128_t *p = &mv->a[i]; uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff; int t; +#ifdef LISA_HASH + t = t_batch[i]; + cr = cr_batch[i]; +#else cr = mm_idx_get(mi, p->x>>8, &t); +#endif if (t == 0) continue; q = &m[k++]; q->q_pos = q_pos, q->q_span = q_span, q->cr = cr, q->n = t, q->seg_id = p->y >> 32; @@ -22,7 +58,16 @@ mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) q->is_tandem = 1; if (i < mv->n - 1 && p->x>>8 == mv->a[i + 1].x>>8) q->is_tandem = 1; } +#ifdef LISA_HASH + free(cr_batch); + free(t_batch); + free(minimizers); + free(lisa_pos); +#endif *n_m_ = k; +#ifdef MANUAL_PROFILING + minimizer_lookup_time += __rdtsc() - lookup_start; +#endif return m; }