diff --git a/Makefile b/Makefile index a49d410..5d685ea 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ CC= gcc #CFLAGS= -g -Wall -Wno-unused-function -mavx2 CFLAGS= -Wall -Wno-unused-function -O2 -mavx2 -DFLAGS= -DSHOW_PERF -DDEBUG_RETURN_VALUE +DFLAGS= -DSHOW_PERF #DFLAGS= -DSHOW_PERF -DDEBUG_OUT -DDEBUG_RETURN_VALUE PROG= sw_perf INCLUDES= diff --git a/ksw_ext_avx2.c b/ksw_ext_avx2.c index c9a1aa0..498b098 100644 --- a/ksw_ext_avx2.c +++ b/ksw_ext_avx2.c @@ -109,13 +109,15 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { fn_vec = _mm256_max_epi16(fn_vec, zero_vec); \ mn_vec = _mm256_max_epi16(mn_vec, zero_vec); \ hn_vec = _mm256_max_epi16(hn_vec, zero_vec); -//int16_t *t_ptr = (int16_t *)&ts_vec; \ - //fprintf(stderr, "D: %d, ibeg: %d, iend: %d, jbeg: %d, jend: %d, %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d \n", \ - // D, ibeg, iend, beg, end, \ - // t_ptr[0], t_ptr[1], t_ptr[2], t_ptr[3], \ - // t_ptr[4], t_ptr[5], t_ptr[6], t_ptr[7], \ - // t_ptr[8], t_ptr[9], t_ptr[10], t_ptr[11], \ - // t_ptr[12], t_ptr[13], t_ptr[14], t_ptr[15]); + +/* int16_t *t_ptr = (int16_t *)&ts_vec; \ + fprintf(stderr, "D: %d, ibeg: %d, iend: %d, jbeg: %d, jend: %d, %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d \n", \ + D, ibeg, iend, beg, end, \ + t_ptr[0], t_ptr[1], t_ptr[2], t_ptr[3], \ + t_ptr[4], t_ptr[5], t_ptr[6], t_ptr[7], \ + t_ptr[8], t_ptr[9], t_ptr[10], t_ptr[11], \ + t_ptr[12], t_ptr[13], t_ptr[14], t_ptr[15]); +*/ // 存储向量化结果 #define SIMD_STORE \ @@ -220,7 +222,7 @@ int ksw_extend_avx2(thread_mem_t *tmem, for (i = 0; i < qlen; ++i) seq[i] = query[qlen - 1 - i]; for (i = 0; i < tlen; ++i) - ref[i + SIMD_WIDTH] = target[tlen - 1 - i]; + ref[i + SIMD_WIDTH - 1] = target[tlen - 1 - i]; } else { @@ -305,16 +307,16 @@ int ksw_extend_avx2(thread_mem_t *tmem, end1 = D; // 闭区间 else end1 = qlen; - // beg1 = MAX(D - window_size, beg1); - // end1 = MIN(D + window_size, end1); + beg1 = MAX(D - window_size, beg1); + end1 = MIN(D + window_size, end1); beg = MAX(beg1, beg); end = MIN(end1, end); - // if (beg > end) - // break; + if (beg > end) + break; - beg = beg1; - end = end1; + // beg = beg1; + // end = end1; iend = D - beg; // ref开始计算的位置，倒序 span = end - beg; @@ -413,7 +415,7 @@ int ksw_extend_avx2(thread_mem_t *tmem, } // 调整计算的边界 - /*for (j = beg; LIKELY(j <= end); ++j) + for (j = beg; LIKELY(j <= end); ++j) { int has_val = hA1[j - 1] | hA2[j]; if (has_val) @@ -425,13 +427,11 @@ int ksw_extend_avx2(thread_mem_t *tmem, int has_val = hA1[j - 1] | hA2[j]; if (has_val) break; - else - hA0[j - 1] = 0; } end = j + 1 <= qlen ? j + 1 : qlen; // beg = 0; // end = qlen; // uncomment this line for debugging - */ + // swap m, h, e, f SWAP_DATA_POINTER; } diff --git a/ksw_ext_avx2_aligned.c b/ksw_ext_avx2_aligned.c index c5bdbe8..e708ea7 100644 --- a/ksw_ext_avx2_aligned.c +++ b/ksw_ext_avx2_aligned.c @@ -85,14 +85,14 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { * h 表示最大值 */ // load向量化数据 -#define SIMD_LOAD \ - __m256i m1 = _mm256_loadu_si256((__m256i *)(&cur_match_arr[j])); \ - __m256i e1 = _mm256_loadu_si256((__m256i *)(&cur_del_arr[j])); \ - __m256i m1j1 = _mm256_loadu_si256((__m256i *)(&cur_match_arr[j - 1])); \ - __m256i f1j1 = _mm256_loadu_si256((__m256i *)(&cur_ins_arr[j - 1])); \ - __m256i h0j1 = _mm256_loadu_si256((__m256i *)(&last_max_arr[j - 1])); \ - __m256i qs_vec = _mm256_loadu_si256((__m256i *)(&read_seq[j])); \ - __m256i ts_vec = _mm256_loadu_si256((__m256i *)(&ref_seq[i])); +#define SIMD_LOAD \ + __m256i m1 = _mm256_loadu_si256((__m256i *)(&mA1[j])); \ + __m256i e1 = _mm256_loadu_si256((__m256i *)(&eA1[j])); \ + __m256i m1j1 = _mm256_loadu_si256((__m256i *)(&mA1[j - 1])); \ + __m256i f1j1 = _mm256_loadu_si256((__m256i *)(&fA1[j - 1])); \ + __m256i h0j1 = _mm256_loadu_si256((__m256i *)(&hA0[j - 1])); \ + __m256i qs_vec = _mm256_loadu_si256((__m256i *)(&seq[j])); \ + __m256i ts_vec = _mm256_loadu_si256((__m256i *)(&ref[i])); // 比对ref和seq的序列，计算罚分 #define SIMD_CMP_SEQ \ @@ -130,74 +130,74 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { mn_vec = _mm256_max_epi16(mn_vec, zero_vec); \ hn_vec = _mm256_max_epi16(hn_vec, zero_vec); -#define SIMD_STORE \ - max_vec = _mm256_max_epu8(max_vec, hn_vec); \ - _mm256_storeu_si256((__m256i *)&next_del_arr[j], en_vec); \ - _mm256_storeu_si256((__m256i *)&next_ins_arr[j], fn_vec); \ - _mm256_storeu_si256((__m256i *)&next_match_arr[j], mn_vec); \ - _mm256_storeu_si256((__m256i *)&next_max_arr[j], hn_vec); +#define SIMD_STORE \ + max_vec = _mm256_max_epu8(max_vec, hn_vec); \ + _mm256_storeu_si256((__m256i *)&eA2[j], en_vec); \ + _mm256_storeu_si256((__m256i *)&fA2[j], fn_vec); \ + _mm256_storeu_si256((__m256i *)&mA2[j], mn_vec); \ + _mm256_storeu_si256((__m256i *)&hA2[j], hn_vec); // 去除多余的部分 -#define SIMD_REMOVE_EXTRA \ - en_vec = _mm256_and_si256(en_vec, h_vec_mask[read_end_pos - j]); \ - fn_vec = _mm256_and_si256(fn_vec, h_vec_mask[read_end_pos - j]); \ - mn_vec = _mm256_and_si256(mn_vec, h_vec_mask[read_end_pos - j]); \ - hn_vec = _mm256_and_si256(hn_vec, h_vec_mask[read_end_pos - j]); +#define SIMD_REMOVE_EXTRA \ + en_vec = _mm256_and_si256(en_vec, h_vec_mask[end - j]); \ + fn_vec = _mm256_and_si256(fn_vec, h_vec_mask[end - j]); \ + mn_vec = _mm256_and_si256(mn_vec, h_vec_mask[end - j]); \ + hn_vec = _mm256_and_si256(hn_vec, h_vec_mask[end - j]); // 找最大值和位置 -#define SIMD_FIND_MAX \ - __m256i cmp_max = _mm256_cmpgt_epi16(max_vec, last_max_vec); \ - uint32_t cmp_result = _mm256_movemask_epi8(cmp_max); \ - if (cmp_result > 0) \ - { \ - max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 2)); \ - max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 4)); \ - max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 6)); \ - max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 8)); \ - max_vec = _mm256_max_epu16(max_vec, _mm256_permute2x128_si256(max_vec, max_vec, 0x01)); \ - int16_t *maxVal = (int16_t *)&max_vec; \ - m = maxVal[0]; \ - for (j = aligned_read_start_pos, i = aligned_ref_end_pos; j <= read_end_pos; j += SIMD_WIDTH, i -= SIMD_WIDTH) \ - { \ - __m256i h2_vec = _mm256_loadu_si256((__m256i *)(&next_max_arr[j])); \ - __m256i vcmp = _mm256_cmpeq_epi16(h2_vec, max_vec); \ - uint32_t mask = _mm256_movemask_epi8(vcmp); \ - if (mask > 0) \ - { \ - int pos = SIMD_WIDTH - 1 - ((__builtin_clz(mask)) >> 1); \ - mj = j - 1 + pos; \ - mi = i - 1 - pos; \ - for (; mj + 1 < qlen && mi + 1 < tlen; mj++, mi++) \ - { \ - if (read_seq[mj + 2] == ref_seq[mi + 1 + SIMD_WIDTH]) \ - { \ - m += base_match_score; \ - } \ - else \ - { \ - break; \ - } \ - } \ - } \ - } \ - last_max_vec = _mm256_set1_epi16(m); \ +#define SIMD_FIND_MAX \ + __m256i cmp_max = _mm256_cmpgt_epi16(max_vec, last_max_vec); \ + uint32_t cmp_result = _mm256_movemask_epi8(cmp_max); \ + if (cmp_result > 0) \ + { \ + max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 2)); \ + max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 4)); \ + max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 6)); \ + max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 8)); \ + max_vec = _mm256_max_epu16(max_vec, _mm256_permute2x128_si256(max_vec, max_vec, 0x01)); \ + int16_t *maxVal = (int16_t *)&max_vec; \ + m = maxVal[0]; \ + for (j = aligned_read_start_pos, i = aligned_ref_end_pos; j <= end; j += SIMD_WIDTH, i -= SIMD_WIDTH) \ + { \ + __m256i h2_vec = _mm256_loadu_si256((__m256i *)(&hA2[j])); \ + __m256i vcmp = _mm256_cmpeq_epi16(h2_vec, max_vec); \ + uint32_t mask = _mm256_movemask_epi8(vcmp); \ + if (mask > 0) \ + { \ + int pos = SIMD_WIDTH - 1 - ((__builtin_clz(mask)) >> 1); \ + mj = j - 1 + pos; \ + mi = i - 1 - pos; \ + for (; mj + 1 < qlen && mi + 1 < tlen; mj++, mi++) \ + { \ + if (seq[mj + 2] == ref[mi + 1 + SIMD_WIDTH]) \ + { \ + m += base_match_score; \ + } \ + else \ + { \ + break; \ + } \ + } \ + } \ + } \ + last_max_vec = _mm256_set1_epi16(m); \ } // 每轮迭代后，交换数组 -#define SWAP_DATA_POINTER \ - int16_t *tmp = last_max_arr; \ - last_max_arr = cur_max_arr; \ - cur_max_arr = next_max_arr; \ - next_max_arr = tmp; \ - tmp = cur_del_arr; \ - cur_del_arr = next_del_arr; \ - next_del_arr = tmp; \ - tmp = cur_ins_arr; \ - cur_ins_arr = next_ins_arr; \ - next_ins_arr = tmp; \ - tmp = cur_match_arr; \ - cur_match_arr = next_match_arr; \ - next_match_arr = tmp; +#define SWAP_DATA_POINTER \ + int16_t *tmp = hA0; \ + hA0 = hA1; \ + hA1 = hA2; \ + hA2 = tmp; \ + tmp = eA1; \ + eA1 = eA2; \ + eA2 = tmp; \ + tmp = fA1; \ + fA1 = fA2; \ + fA2 = tmp; \ + tmp = mA1; \ + mA1 = mA2; \ + mA2 = tmp; int ksw_extend_avx2_aligned(thread_mem_t *tmem, int qlen, // query length 待匹配段碱基的query长度 @@ -220,11 +220,11 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, int *_gscore, // query的端到端匹配得分 int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 { - int16_t *cur_match_arr, *next_match_arr, - *last_max_arr, *cur_max_arr, *next_max_arr, - *cur_del_arr, *next_del_arr, - *cur_ins_arr, *next_ins_arr; // hA0保存上上个col的H，其他的保存上个H E F M - int16_t *read_seq, *ref_seq; + int16_t *mA1, *mA2, + *hA0, *hA1, *hA2, + *eA1, *eA2, + *fA1, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M + int16_t *seq, *ref; uint8_t *mem_addr; int read_size = align_number(qlen * BASE_BYTES + MEM_ALIGN_BYTES); @@ -235,10 +235,10 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, int score_mem_size = score_array_size * TMP_SCORE_ARRAY_NUM; int request_mem_size = read_size + ref_size + score_mem_size + MEM_ALIGN_BYTES * 3; // 左侧内存地址对齐 + 数据向左偏移一个元素 + 末尾SIMD补齐 - int i, ref_start_pos, di, j, read_start_pos, read_end_pos, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; + int i, ibeg, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; int span, beg1, end1; // 边界条件计算 int aligned_read_start_pos, aligned_ref_end_pos; - int ref_end_pos; + int iend; SIMD_INIT; // 初始化simd用的数据 @@ -247,21 +247,21 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, // allocate memory mem_addr = thread_mem_request(tmem, request_mem_size); mem_addr = (void *)align_mem((uint64_t)mem_addr); - ref_seq = (int16_t *)&mem_addr[0]; - read_seq = (int16_t *)(mem_addr + ref_size + SIMD_BYTES - BASE_BYTES); + ref = (int16_t *)&mem_addr[0]; + seq = (int16_t *)(mem_addr + ref_size + SIMD_BYTES - BASE_BYTES); if (extend_left) { for (i = 0; i < qlen; ++i) - read_seq[i + 1] = query[qlen - 1 - i]; + seq[i + 1] = query[qlen - 1 - i]; for (i = 0; i < tlen; ++i) - ref_seq[i + SIMD_WIDTH] = target[tlen - 1 - i]; + ref[i + SIMD_WIDTH - 1] = target[tlen - 1 - i]; } else { for (i = 0; i < qlen; ++i) - read_seq[i + 1] = query[i]; + seq[i + 1] = query[i]; for (i = 0; i < tlen; ++i) - ref_seq[i + SIMD_WIDTH] = target[i]; + ref[i + SIMD_WIDTH - 1] = target[i]; } mem_addr += read_size + ref_size + (SIMD_BYTES - SCORE_BYTES); @@ -270,15 +270,15 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, _mm256_storeu_si256((__m256i *)&mem_addr[i], zero_vec); } - last_max_arr = (int16_t *)&mem_addr[0]; - cur_max_arr = &last_max_arr[score_element_num]; - next_max_arr = &cur_max_arr[score_element_num]; - cur_match_arr = &next_max_arr[score_element_num]; - next_match_arr = &cur_match_arr[score_element_num]; - cur_del_arr = &next_match_arr[score_element_num]; - next_del_arr = &cur_del_arr[score_element_num]; - cur_ins_arr = &next_del_arr[score_element_num]; - next_ins_arr = &cur_ins_arr[score_element_num]; + hA0 = (int16_t *)&mem_addr[0]; + hA1 = &hA0[score_element_num]; + hA2 = &hA1[score_element_num]; + mA1 = &hA2[score_element_num]; + mA2 = &mA1[score_element_num]; + eA1 = &mA2[score_element_num]; + eA2 = &eA1[score_element_num]; + fA1 = &eA2[score_element_num]; + fA2 = &fA1[score_element_num]; // adjust $window_size if it is too large // get the max score @@ -297,10 +297,13 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, max_ie = -1, gscore = -1; ; max_off = 0; - read_start_pos = 1; - read_end_pos = qlen; + beg = 1; + end = qlen; // init init_score - last_max_arr[0] = init_score; // 左上角 + hA0[0] = init_score; // 左上角 + fA1[1] = MAX(0, init_score - (o_ins + e_ins)); + eA2[0] = init_score; + hA1[1] = fA1[1]; if (qlen == 0 || tlen == 0) back_diagnal_num = 0; // 防止意外情况 @@ -310,64 +313,51 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, gscore = 0; } - for (di = 1; LIKELY(di < back_diagnal_num); ++di) + for (D = 1; LIKELY(D < back_diagnal_num); ++D) { - // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 - if (di > tlen) - { - span = MIN(back_diagnal_num - di, window_size); // 计算的窗口，或者说范围 - beg1 = MAX(di - tlen + 1, ((di - window_size) / 2) + 1); - } + if (D < tlen) + beg1 = 1; else - { - span = MIN(di - 1, window_size); - beg1 = MAX(1, ((di - window_size) / 2) + 1); - } - end1 = MIN(qlen, beg1 + span); + beg1 = D - tlen + 1; + if (D < qlen) + end1 = D; // 闭区间 + else + end1 = qlen; + beg1 = MAX(D - window_size, beg1); + end1 = MIN(D + window_size, end1); - // if (read_start_pos < beg1) - // read_start_pos = beg1; - // if (read_end_pos > end1) - // read_end_pos = end1; - // if (read_start_pos > read_end_pos) - // break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug + beg = MAX(beg1, beg); + end = MIN(end1, end); + if (beg > end) + break; - read_start_pos = 1; - read_end_pos = qlen; - ref_end_pos = di - (read_start_pos - 1); // ref开始计算的位置，倒序 - span = read_end_pos - read_start_pos; - ref_start_pos = ref_end_pos - span - 1; // 0开始的ref索引位置 + // beg = beg1; + // end = end1; + + iend = D - beg; // ref开始计算的位置，倒序 + span = end - beg; + ibeg = iend - span; // 0开始的ref索引位置 // 每一轮需要记录的数据 int m = 0, mj = -1, mi = -1; max_vec = zero_vec; - - // 要处理边界 - // 左边界 处理f (insert) - if (ref_start_pos == 0) + // 处理左边界 + if (beg == 1) { - cur_max_arr[read_end_pos] = MAX(0, init_score - (o_ins + e_ins * read_end_pos)); - } - // 上边界 delete - if (read_start_pos == 1) - { - cur_max_arr[0] = MAX(0, init_score - (o_del + e_del * ref_end_pos)); - } - else - { - cur_max_arr[read_start_pos - 1] = 0; - cur_del_arr[read_start_pos - 1] = 0; + hA0[0] = eA2[0]; + mA1[0] = 0; + eA1[0] = MAX(0, init_score - (o_del + e_del * (iend + 1))); } - // aligned_read_start_pos = (read_start_pos >> ALIGN_SHIFT_BITS << ALIGN_SHIFT_BITS) + 1; - // aligned_ref_end_pos = ref_end_pos + (read_start_pos - aligned_read_start_pos); + // aligned_read_start_pos = (beg >> ALIGN_SHIFT_BITS << ALIGN_SHIFT_BITS) + 1; + // aligned_ref_end_pos = iend + (beg - aligned_read_start_pos); - aligned_read_start_pos = read_start_pos; - aligned_ref_end_pos = ref_end_pos; + aligned_read_start_pos = beg; + aligned_ref_end_pos = iend; - // fprintf(stderr, "%d\t%d\n", read_start_pos, aligned_read_start_pos); + // fprintf(stderr, "%d\t%d\n", beg, aligned_read_start_pos); - for (j = aligned_read_start_pos, i = aligned_ref_end_pos; j <= read_end_pos + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH) + for (j = aligned_read_start_pos, i = aligned_ref_end_pos; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH) { // 取数据 SIMD_LOAD; @@ -379,7 +369,7 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, SIMD_STORE; } // 剩下的计算单元 - if (j <= read_end_pos) + if (j <= end) { // 取数据 SIMD_LOAD; @@ -392,16 +382,22 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, // 存储结果 SIMD_STORE; } - + // 处理上边界 + if (ibeg == 0) + { + fA2[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1))); + hA2[end + 1] = fA2[end + 1]; + mA2[end + 1] = 0; + } SIMD_FIND_MAX; // 注意最后跳出循环j的值 - j = read_end_pos + 1; + j = end + 1; if (j == qlen + 1) // 遍历到了query最后一个碱基，此时next_max_arr[qlen]为全局匹配的最大分值 { - max_ie = gscore > next_max_arr[qlen] ? max_ie : ref_start_pos; - gscore = gscore > next_max_arr[qlen] ? gscore : next_max_arr[qlen]; + max_ie = gscore > hA2[qlen] ? max_ie : ibeg; + gscore = gscore > hA2[qlen] ? gscore : hA2[qlen]; } if (m > max) { @@ -410,26 +406,25 @@ int ksw_extend_avx2_aligned(thread_mem_t *tmem, } // 调整计算的边界 - // for (j = read_start_pos; LIKELY(j <= read_end_pos); ++j) - //{ - // int has_val = cur_max_arr[j - 1] | next_max_arr[j]; - // if (has_val) - // { - // break; - // } - //} - // read_start_pos = j; - // - // next_max_arr[read_end_pos + 1] = 0; - // for (j = read_end_pos + 1; LIKELY(j >= read_start_pos); --j) - //{ - // int has_val = cur_max_arr[j - 1] | next_max_arr[j]; - // if (has_val) - // { - // break; - // } - //} - // read_end_pos = j + 1 <= qlen ? j + 1 : qlen; + for (j = beg; LIKELY(j <= end); ++j) + { + int has_val = hA1[j - 1] | hA2[j]; + if (has_val) + { + break; + } + } + beg = j; + + for (j = end + 1; LIKELY(j >= beg); --j) + { + int has_val = hA1[j - 1] | hA2[j]; + if (has_val) + { + break; + } + } + end = j + 1 <= qlen ? j + 1 : qlen; // swap m, h, e, f SWAP_DATA_POINTER; } diff --git a/ksw_ext_avx2_heuristics.c b/ksw_ext_avx2_heuristics.c index aca2f04..9c2050a 100644 --- a/ksw_ext_avx2_heuristics.c +++ b/ksw_ext_avx2_heuristics.c @@ -222,7 +222,7 @@ int ksw_extend_avx2_heuristics(thread_mem_t *tmem, uint8_t *mem; int16_t *qtmem, *vmem; int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH; - int i, iStart, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; + int i, ibeg, iend, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; int Dloop = tlen + qlen; // 循环跳出条件 D从1开始遍历 int span, beg1, end1; // 边界条件计算 int col_size = qlen + 2 + SIMD_WIDTH; @@ -243,14 +243,14 @@ int ksw_extend_avx2_heuristics(thread_mem_t *tmem, for (i = 0; i < qlen; ++i) seq[i] = query[qlen - 1 - i]; for (i = 0; i < tlen; ++i) - ref[i + SIMD_WIDTH] = target[tlen - 1 - i]; + ref[i + SIMD_WIDTH - 1] = target[tlen - 1 - i]; } else { for (i = 0; i < qlen; ++i) seq[i] = query[i]; for (i = 0; i < tlen; ++i) - ref[i + SIMD_WIDTH] = target[i]; + ref[i + SIMD_WIDTH - 1] = target[i]; } vmem = &ref[ref_size]; @@ -288,12 +288,15 @@ int ksw_extend_avx2_heuristics(thread_mem_t *tmem, // DP loop max = init_score, max_i = max_j = -1; max_ie = -1, gscore = -1; - ; + max_off = 0; beg = 1; end = qlen; // init init_score last_max_arr[0] = init_score; // 左上角 + cur_ins_arr[1] = MAX(0, init_score - (o_ins + e_ins)); + next_del_arr[0] = init_score; + cur_max_arr[1] = cur_ins_arr[1]; if (qlen == 0 || tlen == 0) Dloop = 0; // 防止意外情况 @@ -303,56 +306,41 @@ int ksw_extend_avx2_heuristics(thread_mem_t *tmem, gscore = 0; } - int iend; - for (D = 1; LIKELY(D < Dloop); ++D) { - // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 - if (D > tlen) - { - span = MIN(Dloop - D, window_size); // 计算的窗口，或者说范围 - beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1); - } + if (D < tlen) + beg1 = 1; else - { - span = MIN(D - 1, window_size); - beg1 = MAX(1, ((D - window_size) / 2) + 1); - } - end1 = MIN(qlen, beg1 + span); + beg1 = D - tlen + 1; + if (D < qlen) + end1 = D; // 闭区间 + else + end1 = qlen; + beg1 = MAX(D - window_size, beg1); + end1 = MIN(D + window_size, end1); - beg = 1; - end = qlen; + beg = MAX(beg1, beg); + end = MIN(end1, end); + if (beg > end) + break; - // if (beg < beg1) - // beg = beg1; - // if (end > end1) - // end = end1; - // if (beg > end) - // break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug + // beg = beg1; + // end = end1; - iend = D - (beg - 1); // ref开始计算的位置，倒序 + iend = D - beg; // ref开始计算的位置，倒序 span = end - beg; - iStart = iend - span - 1; // 0开始的ref索引位置 + ibeg = iend - span; // 0开始的ref索引位置 // 每一轮需要记录的数据 int m = 0, mj = -1, mi = -1; max_vec = zero_vec; - // 要处理边界 - // 左边界 处理f (insert) - if (iStart == 0) - { - cur_max_arr[end] = MAX(0, init_score - (o_ins + e_ins * end)); - } - // 上边界 + // 处理左边界 if (beg == 1) { - cur_max_arr[0] = MAX(0, init_score - (o_del + e_del * iend)); - } - else - { - cur_max_arr[beg - 1] = 0; - cur_del_arr[beg - 1] = 0; + last_max_arr[0] = next_del_arr[0]; + cur_match_arr[0] = 0; + cur_del_arr[0] = MAX(0, init_score - (o_del + e_del * (iend + 1))); } for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH) @@ -380,6 +368,13 @@ int ksw_extend_avx2_heuristics(thread_mem_t *tmem, // 存储结果 SIMD_STORE; } + // 处理上边界 + if (ibeg == 0) + { + next_ins_arr[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1))); + next_max_arr[end + 1] = next_ins_arr[end + 1]; + next_match_arr[end + 1] = 0; + } SIMD_FIND_MAX; @@ -388,7 +383,7 @@ int ksw_extend_avx2_heuristics(thread_mem_t *tmem, if (j == qlen + 1) { - max_ie = gscore > next_max_arr[qlen] ? max_ie : iStart; + max_ie = gscore > next_max_arr[qlen] ? max_ie : ibeg; gscore = gscore > next_max_arr[qlen] ? gscore : next_max_arr[qlen]; } if (m > max) @@ -397,25 +392,24 @@ int ksw_extend_avx2_heuristics(thread_mem_t *tmem, max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi); } - // for (j = beg; LIKELY(j <= end); ++j) - //{ - // int has_val = cur_max_arr[j - 1] | next_max_arr[j]; - // if (has_val) - // { - // break; - // } - // } - // beg = j; - // next_max_arr[end + 1] = 0; - // for (j = end + 1; LIKELY(j >= beg); --j) - //{ - // int has_val = cur_max_arr[j - 1] | next_max_arr[j]; - // if (has_val) - // { - // break; - // } - // } - // end = j + 1 <= qlen ? j + 1 : qlen; + for (j = beg; LIKELY(j <= end); ++j) + { + int has_val = cur_max_arr[j - 1] | next_max_arr[j]; + if (has_val) + { + break; + } + } + beg = j; + for (j = end + 1; LIKELY(j >= beg); --j) + { + int has_val = cur_max_arr[j - 1] | next_max_arr[j]; + if (has_val) + { + break; + } + } + end = j + 1 <= qlen ? j + 1 : qlen; // swap m, h, e, f SWAP_DATA_POINTER; diff --git a/ksw_ext_avx2_u8.c b/ksw_ext_avx2_u8.c index 8b38cb3..5053051 100644 --- a/ksw_ext_avx2_u8.c +++ b/ksw_ext_avx2_u8.c @@ -243,7 +243,7 @@ int ksw_extend_avx2_u8(thread_mem_t *tmem, for (i = 0; i < qlen; ++i) seq[i] = query[qlen - 1 - i]; for (i = 0; i < tlen; ++i) - ref[i + SIMD_WIDTH] = target[tlen - 1 - i]; + ref[i + SIMD_WIDTH - 1] = target[tlen - 1 - i]; } else { @@ -317,16 +317,16 @@ int ksw_extend_avx2_u8(thread_mem_t *tmem, end1 = D; // 闭区间 else end1 = qlen; - // beg1 = MAX(D - window_size, beg1); - // end1 = MIN(D + window_size, end1); + beg1 = MAX(D - window_size, beg1); + end1 = MIN(D + window_size, end1); beg = MAX(beg1, beg); end = MIN(end1, end); - // if (beg > end) - // break; + if (beg > end) + break; - beg = beg1; - end = end1; + // beg = beg1; + // end = end1; iend = D - beg; // ref开始计算的位置，倒序 span = end - beg; @@ -394,22 +394,22 @@ int ksw_extend_avx2_u8(thread_mem_t *tmem, } // 调整计算的边界 - // for (j = beg; LIKELY(j <= end); ++j) - //{ - // int has_val = hA1[j - 1] | hA2[j]; - // if (has_val) - // break; - //} - // beg = j; - // for (j = end + 1; LIKELY(j >= beg); --j) - //{ - // int has_val = hA1[j - 1] | hA2[j]; - // if (has_val) - // break; - // else - // hA0[j - 1] = 0; - //} - // end = j + 1 <= qlen ? j + 1 : qlen; + for (j = beg; LIKELY(j <= end); ++j) + { + int has_val = hA1[j - 1] | hA2[j]; + if (has_val) + break; + } + beg = j; + for (j = end + 1; LIKELY(j >= beg); --j) + { + int has_val = hA1[j - 1] | hA2[j]; + if (has_val) + break; + else + hA0[j - 1] = 0; + } + end = j + 1 <= qlen ? j + 1 : qlen; // swap m, h, e, f SWAP_DATA_POINTER; diff --git a/ksw_ext_avx2_u8_aligned.c b/ksw_ext_avx2_u8_aligned.c index ba4cdc4..1aad002 100644 --- a/ksw_ext_avx2_u8_aligned.c +++ b/ksw_ext_avx2_u8_aligned.c @@ -250,7 +250,7 @@ int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, { uint8_t *mA1, *mA2, *hA0, *hA1, *eA1, *fA1, *hA2, *eA2, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M uint8_t *read_seq, *ref_seq; - int i, iStart, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; + int i, ibeg, iend, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; int span, beg1, end1; // 边界条件计算 uint8_t *mem_addr; @@ -276,14 +276,14 @@ int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, for (i = 0; i < qlen; ++i) read_seq[i + 1] = query[qlen - 1 - i]; for (i = 0; i < tlen; ++i) - ref_seq[i + SIMD_WIDTH] = target[tlen - 1 - i]; + ref_seq[i + SIMD_WIDTH - 1] = target[tlen - 1 - i]; } else { for (i = 0; i < qlen; ++i) read_seq[i + 1] = query[i]; for (i = 0; i < tlen; ++i) - ref_seq[i + SIMD_WIDTH] = target[i]; + ref_seq[i + SIMD_WIDTH - 1] = target[i]; } mem_addr += read_size + ref_size; @@ -318,12 +318,15 @@ int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, // DP loop max = init_score, max_i = max_j = -1; max_ie = -1, gscore = -1; - ; + max_off = 0; beg = 1; end = qlen; // init init_score hA0[0] = init_score; // 左上角 + fA1[1] = MAX(0, init_score - (o_ins + e_ins)); + eA2[0] = init_score; + hA1[1] = fA1[1]; if (qlen == 0 || tlen == 0) back_diagnal_num = 0; // 防止意外情况 @@ -333,55 +336,42 @@ int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, gscore = 0; } - int iend; - for (D = 1; LIKELY(D < back_diagnal_num); ++D) { // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 - if (D > tlen) - { - span = MIN(back_diagnal_num - D, window_size); - beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1); - } + if (D < tlen) + beg1 = 1; else - { - span = MIN(D - 1, window_size); - beg1 = MAX(1, ((D - window_size) / 2) + 1); - } - end1 = MIN(qlen, beg1 + span); + beg1 = D - tlen + 1; + if (D < qlen) + end1 = D; // 闭区间 + else + end1 = qlen; + beg1 = MAX(D - window_size, beg1); + end1 = MIN(D + window_size, end1); - // if (beg < beg1) - // beg = beg1; - // if (end > end1) - // end = end1; - // if (beg > end) - // break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug + beg = MAX(beg1, beg); + end = MIN(end1, end); + if (beg > end) + break; - beg = 1; - end = qlen; - iend = D - (beg - 1); // ref开始计算的位置，倒序 + // beg = beg1; + // end = end1; + + iend = D - beg; // ref开始计算的位置，倒序 span = end - beg; - iStart = iend - span - 1; // 0开始的ref索引位置 + ibeg = iend - span; // 0开始的ref索引位置 // 每一轮需要记录的数据 int m = 0, mj = -1, mi = -1; max_vec = zero_vec; - // 要处理边界 - // 左边界 处理f (insert) - if (iStart == 0) - { - hA1[end] = MAX(0, init_score - (o_ins + e_ins * end)); - } - // 上边界 + // 处理左边界 if (beg == 1) { - hA1[0] = MAX(0, init_score - (o_del + e_del * iend)); - } - else - { - hA1[beg - 1] = 0; - eA1[beg - 1] = 0; + hA0[0] = eA2[0]; + mA1[0] = 0; + eA1[0] = MAX(0, init_score - (o_del + e_del * (iend + 1))); } for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH) @@ -410,6 +400,13 @@ int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, SIMD_STORE; } + // 处理上边界 + if (ibeg == 0) + { + fA2[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1))); + hA2[end + 1] = fA2[end + 1]; + mA2[end + 1] = 0; + } SIMD_FIND_MAX; // 注意最后跳出循环j的值 @@ -417,7 +414,7 @@ int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, if (j == qlen + 1) { - max_ie = gscore > hA2[qlen] ? max_ie : iStart; + max_ie = gscore > hA2[qlen] ? max_ie : ibeg; gscore = gscore > hA2[qlen] ? gscore : hA2[qlen]; } @@ -428,21 +425,20 @@ int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, } // 调整计算的边界 - // for (j = beg; LIKELY(j <= end); ++j) - //{ - // int has_val = hA1[j - 1] | hA2[j]; - // if (has_val) - // break; - //} - // beg = j; - // hA2[end + 1] = 0; - // for (j = end + 1; LIKELY(j >= beg); --j) - //{ - // int has_val = hA1[j - 1] | hA2[j]; - // if (has_val) - // break; - //} - // end = j + 1 <= qlen ? j + 1 : qlen; + for (j = beg; LIKELY(j <= end); ++j) + { + int has_val = hA1[j - 1] | hA2[j]; + if (has_val) + break; + } + beg = j; + for (j = end + 1; LIKELY(j >= beg); --j) + { + int has_val = hA1[j - 1] | hA2[j]; + if (has_val) + break; + } + end = j + 1 <= qlen ? j + 1 : qlen; // swap m, h, e, f SWAP_DATA_POINTER; diff --git a/ksw_ext_avx2_u8_heuristics.c b/ksw_ext_avx2_u8_heuristics.c index 72d4690..4a85d3d 100644 --- a/ksw_ext_avx2_u8_heuristics.c +++ b/ksw_ext_avx2_u8_heuristics.c @@ -243,7 +243,7 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, uint8_t *seq, *ref; uint8_t *mem, *qtmem, *vmem; int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH; - int i, iStart, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; + int i, ibeg, iend, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; int Dloop = tlen + qlen; // 循环跳出条件 int span, beg1, end1; // 边界条件计算 int col_size = qlen + 2 + SIMD_WIDTH; @@ -265,14 +265,14 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, for (i = 0; i < qlen; ++i) seq[i] = query[qlen - 1 - i]; for (i = 0; i < tlen; ++i) - ref[i + SIMD_WIDTH] = target[tlen - 1 - i]; + ref[i + SIMD_WIDTH - 1] = target[tlen - 1 - i]; } else { for (i = 0; i < qlen; ++i) seq[i] = query[i]; for (i = 0; i < tlen; ++i) - ref[i + SIMD_WIDTH] = target[i]; + ref[i + SIMD_WIDTH - 1] = target[i]; } vmem = &ref[ref_size]; @@ -311,12 +311,15 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, // DP loop max = init_score, max_i = max_j = -1; max_ie = -1, gscore = -1; - ; + max_off = 0; beg = 1; end = qlen; // init init_score hA0[0] = init_score; // 左上角 + fA1[1] = MAX(0, init_score - (o_ins + e_ins)); + eA2[0] = init_score; + hA1[1] = fA1[1]; if (qlen == 0 || tlen == 0) Dloop = 0; // 防止意外情况 @@ -325,62 +328,42 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, max_ie = 0; gscore = 0; } - - int iend; -#ifdef SHOW_PERF -// time_bsw_init += get_mseconds() - start_time; -#endif - for (D = 1; LIKELY(D < Dloop); ++D) { -#ifdef SHOW_PERF -// start_time = get_mseconds(); -#endif // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 - if (D > tlen) - { - span = MIN(Dloop - D, window_size); - beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1); - } + if (D < tlen) + beg1 = 1; else - { - span = MIN(D - 1, window_size); - beg1 = MAX(1, ((D - window_size) / 2) + 1); - } - end1 = MIN(qlen, beg1 + span); + beg1 = D - tlen + 1; + if (D < qlen) + end1 = D; // 闭区间 + else + end1 = qlen; + beg1 = MAX(D - window_size, beg1); + end1 = MIN(D + window_size, end1); - // if (beg < beg1) - // beg = beg1; - // if (end > end1) - // end = end1; - // if (beg > end) - // break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug + beg = MAX(beg1, beg); + end = MIN(end1, end); + if (beg > end) + break; - beg = 1; - end = qlen; - iend = D - (beg - 1); // ref开始计算的位置，倒序 + // beg = beg1; + // end = end1; + + iend = D - beg; // ref开始计算的位置，倒序 span = end - beg; - iStart = iend - span - 1; // 0开始的ref索引位置 + ibeg = iend - span; // 0开始的ref索引位置 // 每一轮需要记录的数据 int m = 0, mj = -1, mi = -1; max_vec = zero_vec; - // 要处理边界 - // 左边界 处理f (insert) - if (iStart == 0) - { - hA1[end] = MAX(0, init_score - (o_ins + e_ins * end)); - } - // 上边界 + // 处理左边界 if (beg == 1) { - hA1[0] = MAX(0, init_score - (o_del + e_del * iend)); - } - else - { - hA1[beg - 1] = 0; - eA1[beg - 1] = 0; + hA0[0] = eA2[0]; + mA1[0] = 0; + eA1[0] = MAX(0, init_score - (o_del + e_del * (iend + 1))); } for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH) @@ -408,6 +391,13 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, // 存储结果 SIMD_STORE; } + // 处理上边界 + if (ibeg == 0) + { + fA2[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1))); + hA2[end + 1] = fA2[end + 1]; + mA2[end + 1] = 0; + } SIMD_FIND_MAX; @@ -416,7 +406,7 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, if (j == qlen + 1) { - max_ie = gscore > hA2[qlen] ? max_ie : iStart; + max_ie = gscore > hA2[qlen] ? max_ie : ibeg; gscore = gscore > hA2[qlen] ? gscore : hA2[qlen]; } @@ -427,21 +417,20 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, } // 调整计算的边界 - // for (j = beg; LIKELY(j <= end); ++j) - //{ - // int has_val = hA1[j - 1] | hA2[j]; - // if (has_val) - // break; - //} - // beg = j; - // hA2[end + 1] = 0; - // for (j = end + 1; LIKELY(j >= beg); --j) - //{ - // int has_val = hA1[j - 1] | hA2[j]; - // if (has_val) - // break; - //} - // end = j + 1 <= qlen ? j + 1 : qlen; + for (j = beg; LIKELY(j <= end); ++j) + { + int has_val = hA1[j - 1] | hA2[j]; + if (has_val) + break; + } + beg = j; + for (j = end + 1; LIKELY(j >= beg); --j) + { + int has_val = hA1[j - 1] | hA2[j]; + if (has_val) + break; + } + end = j + 1 <= qlen ? j + 1 : qlen; // swap m, h, e, f SWAP_DATA_POINTER; diff --git a/ksw_ext_normal.c b/ksw_ext_normal.c index cc404c7..3492ce5 100644 --- a/ksw_ext_normal.c +++ b/ksw_ext_normal.c @@ -25,11 +25,11 @@ int ksw_extend_normal(thread_mem_t *tmem, int qlen, const uint8_t *query, int tl int8_t *qp; // query profile int i, j, k, oe_del = o_del + e_del, oe_ins = o_ins + e_ins, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; assert(h0 > 0); - // qp = malloc(qlen * m); - // eh = calloc(qlen + 1, 8); - qp = thread_mem_request(tmem, qlen * m); - eh = thread_mem_request_and_clean(tmem, (qlen + 1) * 8); - // generate the query profile + qp = malloc(qlen * m); + eh = calloc(qlen + 1, 8); + // qp = thread_mem_request(tmem, qlen * m); + // eh = thread_mem_request_and_clean(tmem, (qlen + 1) * 9); + // generate the query profile for (k = i = 0; k < m; ++k) { const int8_t *p = &mat[k * m]; @@ -77,14 +77,14 @@ int ksw_extend_normal(thread_mem_t *tmem, int qlen, const uint8_t *query, int tl // 对于target第i个字符，query中每个字符的分值，只有匹配和不匹配 int8_t *q = &qp[target[i] * qlen]; // apply the band and the constraint (if provided) - // if (beg < i - w) // 检查开始点是否可以缩小一些 - // beg = i - w; - // if (end > i + w + 1) // 检查终点是否可以缩小，使得整体的遍历范围缩小 - // end = i + w + 1; - // if (end > qlen) // 终点不超过query长度 - // end = qlen; - beg = 0; - end = qlen; + if (beg < i - w) // 检查开始点是否可以缩小一些 + beg = i - w; + if (end > i + w + 1) // 检查终点是否可以缩小，使得整体的遍历范围缩小 + end = i + w + 1; + if (end > qlen) // 终点不超过query长度 + end = qlen; + // beg = 0; + // end = qlen; // compute the first column if (beg == 0) { @@ -142,8 +142,8 @@ int ksw_extend_normal(thread_mem_t *tmem, int qlen, const uint8_t *query, int tl max_ie = gscore > h1 ? max_ie : i; // max_ie表示取得全局最大分值时，target字符串的位置 gscore = gscore > h1 ? gscore : h1; } - // if (m == 0) // 遍历完query之后，当前轮次的最大分值为0，则跳出循环 - // break; + if (m == 0) // 遍历完query之后，当前轮次的最大分值为0，则跳出循环 + break; if (m > max) // 当前轮最大分值大于之前的最大分值 { max = m, max_i = i, max_j = mj; // 更新取得最大值的target和query的位置 @@ -163,19 +163,19 @@ int ksw_extend_normal(thread_mem_t *tmem, int qlen, const uint8_t *query, int tl } } // update beg and end for the next round - // for (j = beg; LIKELY(j < end) && eh[j].h == 0 && eh[j].e == 0; ++j) - // ; - // beg = j; - // for (j = end; LIKELY(j >= beg) && eh[j].h == 0 && eh[j].e == 0; --j) - // ; - // end = j + 2 < qlen ? j + 2 : qlen; // 剪枝没考虑f，即insert - // beg = 0, end = qlen; // uncomment this line for debugging - // fprintf(stderr, "\n"); - // fprintf(stderr, "%d\n", end); + for (j = beg; LIKELY(j < end) && eh[j].h == 0 && eh[j].e == 0; ++j) + ; + beg = j; + for (j = end; LIKELY(j >= beg) && eh[j].h == 0 && eh[j].e == 0; --j) + ; + end = j + 2 < qlen ? j + 2 : qlen; // 剪枝没考虑f，即insert + // beg = 0, end = qlen; // uncomment this line for debugging + // fprintf(stderr, "\n"); + // fprintf(stderr, "%d\n", end); } - // free(eh); - // free(qp); - thread_mem_release(tmem, qlen * m + (qlen + 1) * 8); + free(eh); + free(qp); + // thread_mem_release(tmem, qlen * m + (qlen + 1) * 9); if (_qle) *_qle = max_j + 1; if (_tle) diff --git a/main.c b/main.c index 137f592..f5adbcc 100644 --- a/main.c +++ b/main.c @@ -242,17 +242,17 @@ int main(int argc, char *argv[]) PERFORMANCE_TEST_NORMAL(0, ksw_extend_normal); // avx2 - PERFORMANCE_TEST_AVX2(1, ksw_extend_avx2); + // PERFORMANCE_TEST_AVX2(1, ksw_extend_avx2); // avx2 heuristics - // PERFORMANCE_TEST_AVX2(2, ksw_extend_avx2_heuristics); - // // avx2 mem aligned - // PERFORMANCE_TEST_AVX2(3, ksw_extend_avx2_aligned); - // - // // avx2 u8 - PERFORMANCE_TEST_AVX2(4, ksw_extend_avx2_u8); - // // avx2 u8 heuristics - // PERFORMANCE_TEST_AVX2(5, ksw_extend_avx2_u8_heuristics); - // // avx2 u8 mem aligned + // PERFORMANCE_TEST_AVX2(2, ksw_extend_avx2_heuristics); + // avx2 mem aligned + // PERFORMANCE_TEST_AVX2(3, ksw_extend_avx2_aligned); + + // avx2 u8 + // PERFORMANCE_TEST_AVX2(4, ksw_extend_avx2_u8); + // avx2 u8 heuristics + // PERFORMANCE_TEST_AVX2(5, ksw_extend_avx2_u8_heuristics); + // avx2 u8 mem aligned // PERFORMANCE_TEST_AVX2(6, ksw_extend_avx2_u8_aligned); } diff --git a/thread_mem.c b/thread_mem.c index 3a5a930..38722f0 100644 --- a/thread_mem.c +++ b/thread_mem.c @@ -38,6 +38,7 @@ void thread_mem_init_alloc(thread_mem_t *tmem, size_t byte_cnt) // 请求内存 void *thread_mem_request(thread_mem_t *tmem, size_t byte_cnt) { + // fprintf(stderr, "capacity:%ld, occupied: %ld, byte_cnt: %ld\n", tmem->capacity, tmem->occupied, byte_cnt); void *ret_mem = 0; if (tmem == 0) {