commit 2904e87dee63faef8e05c75215a3456175515840 Author: zzh Date: Thu Aug 10 15:28:45 2023 +0800 初始化仓库，已经实现了normal和avx2的sw，并进行了性能测试 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..9f84111 --- /dev/null +++ b/.gitignore @@ -0,0 +1,60 @@ +*.[oa] +sw_perf +test +test64 +.*.swp +Makefile.bak +bwamem-lite +# ---> C +# Prerequisites +*.d + +# Object files +*.o +*.ko +*.obj +*.elf + +# Linker output +*.ilk +*.map +*.exp + +# Precompiled Headers +*.gch +*.pch + +# Libraries +*.lib +*.a +*.la +*.lo + +# Shared objects (inc. Windows DLLs) +*.dll +*.so +*.so.* +*.dylib + +# Executables +*.exe +*.out +*.app +*.i*86 +*.x86_64 +*.hex + +# Debug files +*.dSYM/ +*.su +*.idb +*.pdb + +# Kernel Module Compile Results +*.mod* +*.cmd +.tmp_versions/ +modules.order +Module.symvers +Mkfile.old +dkms.conf diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..5e68af2 --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,19 @@ +{ + // 使用 IntelliSense 了解相关属性。 + // 悬停以查看现有属性的描述。 + // 欲了解更多信息，请访问: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "sw-perf", + "preLaunchTask": "Build", + "type": "cppdbg", + "request": "launch", + "program": "${workspaceRoot}/sw_perf", + "args": [ + "all" + ], + "cwd": "${workspaceFolder}", // 当前工作路径：当前文件所在的工作空间 + } + ] +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..b03a764 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "files.associations": { + "functional": "c" + } +} \ No newline at end of file diff --git a/.vscode/tasks.json b/.vscode/tasks.json new file mode 100644 index 0000000..f76ae19 --- /dev/null +++ b/.vscode/tasks.json @@ -0,0 +1,17 @@ +{ + // See https://go.microsoft.com/fwlink/?LinkId=733558 + // for the documentation about the tasks.json format + "version": "2.0.0", + "tasks": [ + { + "label": "Build", + "type": "shell", + "command": "make clean; make -j 16", + "problemMatcher": [], + "group": { + "kind": "build", + "isDefault": true + } + } + ] +} \ No newline at end of file diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..37c830b --- /dev/null +++ b/Makefile @@ -0,0 +1,31 @@ +CC= gcc +#CFLAGS= -g -Wall -Wno-unused-function +CFLAGS= -Wall -Wno-unused-function -O2 -mavx2 +DFLAGS= -DSHOW_PERF +OBJS= ksw_normal.o ksw_avx2.o ksw_cuda.o ksw_avx2_u8.o +PROG= sw_perf +INCLUDES= +LIBS= +SUBDIRS= . + +ifeq ($(shell uname -s),Linux) + LIBS += -lrt +endif + +.SUFFIXES:.c .o .cc + +.c.o: + $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) $(CPPFLAGS) $< -o $@ + +all:$(PROG) + +sw_perf:$(OBJS) main.o + $(CC) $(CFLAGS) $(LDFLAGS) $(OBJS) main.o -o $@ -L. $(LIBS) + +clean: + rm -f *.o a.out $(PROG) *~ *.a + +depend: + ( LC_ALL=C ; export LC_ALL; makedepend -Y -- $(CFLAGS) $(DFLAGS) $(CPPFLAGS) -- *.c ) + +# DO NOT DELETE THIS LINE -- make depend depends on it. \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..e69de29 diff --git a/ksw_avx2.c b/ksw_avx2.c new file mode 100644 index 0000000..233d672 --- /dev/null +++ b/ksw_avx2.c @@ -0,0 +1,612 @@ +#include +#include +#include +#include +#include +#include +#include + +#ifdef __GNUC__ +#define LIKELY(x) __builtin_expect((x), 1) +#define UNLIKELY(x) __builtin_expect((x), 0) +#else +#define LIKELY(x) (x) +#define UNLIKELY(x) (x) +#endif + +#undef MAX +#undef MIN +#define MAX(x, y) ((x) > (y) ? (x) : (y)) +#define MIN(x, y) ((x) < (y) ? (x) : (y)) +#define SIMD_WIDTH 16 + +int ksw_extend2_origin(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int is_left, int m, const int8_t *mat, int o_del, int e_del, + int o_ins, int e_ins, int w, int end_bonus, int zdrop, int h0, int *_qle, int *_tle, int *_gtle, int *_gscore, int *_max_off); + +static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { + {0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0}, + {0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}}; + +// const int permute_mask = _MM_SHUFFLE(0, 1, 2, 3); +#define permute_mask _MM_SHUFFLE(0, 1, 2, 3) +// 初始化变量 +#define SIMD_INIT \ + int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \ + __m256i zero_vec; \ + __m256i max_vec; \ + __m256i oe_del_vec; \ + __m256i oe_ins_vec; \ + __m256i e_del_vec; \ + __m256i e_ins_vec; \ + __m256i h_vec_mask[SIMD_WIDTH]; \ + zero_vec = _mm256_setzero_si256(); \ + oe_del_vec = _mm256_set1_epi16(-oe_del); \ + oe_ins_vec = _mm256_set1_epi16(-oe_ins); \ + e_del_vec = _mm256_set1_epi16(-e_del); \ + e_ins_vec = _mm256_set1_epi16(-e_ins); \ + __m256i match_sc_vec = _mm256_set1_epi16(a); \ + __m256i mis_sc_vec = _mm256_set1_epi16(-b); \ + __m256i amb_sc_vec = _mm256_set1_epi16(-1); \ + __m256i amb_vec = _mm256_set1_epi16(4); \ + for (i = 0; i < SIMD_WIDTH; ++i) \ + h_vec_mask[i] = _mm256_loadu_si256((__m256i *)(&h_vec_int_mask[i])); + +/* + * e 表示当前ref的碱基被删除 + * f 表示当前seq的碱基插入 + * m 表示当前碱基匹配（可以相等，也可以不想等） + * h 表示最大值 + */ +// load向量化数据 +#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 - 1])); \ + __m256i ts_vec = _mm256_loadu_si256((__m256i *)(&ref[i])); + +// 比对ref和seq的序列，计算罚分 +#define SIMD_CMP_SEQ \ + ts_vec = _mm256_permute4x64_epi64(ts_vec, permute_mask); \ + ts_vec = _mm256_shufflelo_epi16(ts_vec, permute_mask); \ + ts_vec = _mm256_shufflehi_epi16(ts_vec, permute_mask); \ + __m256i match_mask_vec = _mm256_cmpeq_epi16(qs_vec, ts_vec); \ + __m256i mis_score_vec = _mm256_andnot_si256(match_mask_vec, mis_sc_vec); \ + __m256i score_vec = _mm256_and_si256(match_sc_vec, match_mask_vec); \ + score_vec = _mm256_or_si256(score_vec, mis_score_vec); \ + __m256i q_amb_mask_vec = _mm256_cmpeq_epi16(qs_vec, amb_vec); \ + __m256i t_amb_mask_vec = _mm256_cmpeq_epi16(ts_vec, amb_vec); \ + __m256i amb_mask_vec = _mm256_or_si256(q_amb_mask_vec, t_amb_mask_vec); \ + score_vec = _mm256_andnot_si256(amb_mask_vec, score_vec); \ + __m256i amb_score_vec = _mm256_and_si256(amb_mask_vec, amb_sc_vec); \ + score_vec = _mm256_or_si256(score_vec, amb_score_vec); + +// 向量化计算h, e, f, m +#define SIMD_COMPUTE \ + __m256i en_vec0 = _mm256_add_epi16(m1, oe_del_vec); \ + __m256i en_vec1 = _mm256_add_epi16(e1, e_del_vec); \ + __m256i en_vec = _mm256_max_epi16(en_vec0, en_vec1); \ + __m256i fn_vec0 = _mm256_add_epi16(m1j1, oe_ins_vec); \ + __m256i fn_vec1 = _mm256_add_epi16(f1j1, e_ins_vec); \ + __m256i fn_vec = _mm256_max_epi16(fn_vec0, fn_vec1); \ + __m256i mn_vec0 = _mm256_add_epi16(h0j1, score_vec); \ + __m256i mn_mask = _mm256_cmpgt_epi16(h0j1, zero_vec); \ + __m256i mn_vec = _mm256_and_si256(mn_vec0, mn_mask); \ + __m256i hn_vec0 = _mm256_max_epi16(en_vec, fn_vec); \ + __m256i hn_vec = _mm256_max_epi16(hn_vec0, mn_vec); \ + en_vec = _mm256_max_epi16(en_vec, zero_vec); \ + 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); + +// 存储向量化结果 +#define SIMD_STORE \ + max_vec = _mm256_max_epi16(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[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 \ + 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]; \ + if (m > 0) \ + { \ + for (j = beg, i = iend; 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; \ + } \ + } \ + } + +// 每轮迭代后，交换数组 +#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_avx2(int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int is_left, // 是不是向左扩展 + int m, // 碱基种类 (5) + const int8_t *mat, // 每个位置的query和target的匹配得分 m*m + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int a, // 碱基match时的分数 + int b, // 碱基mismatch时的惩罚分数（正数） + int w, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, + int zdrop, + int h0, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 +{ + // ksw_extend2_origin + // return ksw_extend2_origin(qlen, query, tlen, target, is_left, m, mat, o_del, e_del, o_ins, e_ins, w, end_bonus, zdrop, h0, _qle, _tle, _gtle, _gscore, _max_off); + // fprintf(stderr, "qlen: %d, tlen: %d\n", qlen, tlen); + + // if (qlen * a + h0 < 255) + // return ksw_extend2_avx2_u8(qlen, query, tlen, target, is_left, m, mat, o_del, e_del, o_ins, e_ins, a, b, w, end_bonus, zdrop, h0, _qle, _tle, _gtle, _gscore, _max_off); + + int16_t *mA, *hA, *eA, *fA, *mA1, *mA2, *hA0, *hA1, *eA1, *fA1, *hA2, *eA2, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M + int16_t *seq, *ref; + uint8_t *mem; + int16_t *qtmem, *vmem; + int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH; + int i, iStart, D, j, k, 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; + int val_mem_size = (col_size * 9 * 2 + 31) >> 5 << 5; // 32字节的整数倍 + int mem_size = (seq_size + ref_size) * 2 + val_mem_size; + + SIMD_INIT; // 初始化simd用的数据 + + assert(h0 > 0); + + // allocate memory + mem = malloc(mem_size); + qtmem = (int16_t *)&mem[0]; + seq = &qtmem[0]; + ref = &qtmem[seq_size]; + if (is_left) + { + 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]; + } + else + { + for (i = 0; i < qlen; ++i) + seq[i] = query[i]; + for (i = 0; i < tlen; ++i) + ref[i + SIMD_WIDTH] = target[i]; + } + + vmem = &ref[ref_size]; + for (i = 0; i < (val_mem_size >> 1); i += SIMD_WIDTH) + { + _mm256_storeu_si256((__m256i *)&vmem[i], zero_vec); + } + hA = &vmem[0]; + mA = &vmem[col_size * 3]; + eA = &vmem[col_size * 5]; + fA = &vmem[col_size * 7]; + + hA0 = &hA[0]; + hA1 = &hA[col_size]; + hA2 = &hA1[col_size]; + mA1 = &mA[0]; + mA2 = &mA[col_size]; + eA1 = &eA[0]; + eA2 = &eA[col_size]; + fA1 = &fA[0]; + fA2 = &fA[col_size]; + + // adjust $w if it is too large + k = m * m; + // get the max score + for (i = 0, max = 0; i < k; ++i) + max = max > mat[i] ? max : mat[i]; + max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.); + max_ins = max_ins > 1 ? max_ins : 1; + w = w < max_ins ? w : max_ins; + max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.); + max_del = max_del > 1 ? max_del : 1; + w = w < max_del ? w : max_del; // TODO: is this necessary? + if (tlen < qlen) + w = MIN(tlen - 1, w); + + // DP loop + max = h0, max_i = max_j = -1; + max_ie = -1, gscore = -1; + ; + max_off = 0; + beg = 1; + end = qlen; + // init h0 + hA0[0] = h0; // 左上角 + + if (qlen == 0 || tlen == 0) + Dloop = 0; // 防止意外情况 + if (w >= qlen) + { + max_ie = 0; + gscore = 0; + } + + int m_last = 0; + int iend; + + for (D = 1; LIKELY(D < Dloop); ++D) + { + // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 + if (D > tlen) + { + span = MIN(Dloop - D, w); + beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1); + } + else + { + span = MIN(D - 1, w); + beg1 = MAX(1, ((D - w) / 2) + 1); + } + end1 = MIN(qlen, beg1 + span); + + if (beg < beg1) + beg = beg1; + if (end > end1) + end = end1; + if (beg > end) + break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug + + iend = D - (beg - 1); // ref开始计算的位置，倒序 + span = end - beg; + iStart = iend - span - 1; // 0开始的ref索引位置 + + // 每一轮需要记录的数据 + int m = 0, mj = -1, mi = -1; + max_vec = zero_vec; + + // 要处理边界 + // 左边界 处理f (insert) + if (iStart == 0) + { + hA1[end] = MAX(0, h0 - (o_ins + e_ins * end)); + } + // 上边界 + if (beg == 1) + { + hA1[0] = MAX(0, h0 - (o_del + e_del * iend)); + } + else + { + hA1[beg - 1] = 0; + eA1[beg - 1] = 0; + } + + for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH) + { + // 取数据 + SIMD_LOAD; + // 比对seq，计算罚分 + SIMD_CMP_SEQ; + // 计算 + SIMD_COMPUTE; + // 存储结果 + SIMD_STORE; + } + // 剩下的计算单元 + if (j <= end) + { + // 取数据 + SIMD_LOAD; + // 比对seq，计算罚分 + SIMD_CMP_SEQ; + // 计算 + SIMD_COMPUTE; + // 去除多余计算的部分 + SIMD_REMOVE_EXTRA; + // 存储结果 + SIMD_STORE; + } + + SIMD_FIND_MAX; + + // 注意最后跳出循环j的值 + j = end + 1; + + if (j == qlen + 1) + { + max_ie = gscore > hA2[qlen] ? max_ie : iStart; + gscore = gscore > hA2[qlen] ? gscore : hA2[qlen]; + } + if (m == 0 && m_last == 0) + break; // 一定要注意，斜对角遍历和按列遍历的不同点 + if (m > max) + { + max = m, max_i = mi, max_j = mj; + max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi); + } + else if (zdrop > 0) + { + if (mi - max_i > mj - max_j) + { + if (max - m - ((mi - max_i) - (mj - max_j)) * e_del > zdrop) + break; + } + else + { + if (max - m - ((mj - max_j) - (mi - max_i)) * e_ins > zdrop) + break; + } + } + + // 调整计算的边界 + 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; + + m_last = m; + // swap m, h, e, f + SWAP_DATA_POINTER; + } + + free(mem); + if (_qle) + *_qle = max_j + 1; + if (_tle) + *_tle = max_i + 1; + if (_gtle) + *_gtle = max_ie + 1; + if (_gscore) + *_gscore = gscore; + if (_max_off) + *_max_off = max_off; + return max; +} + +typedef struct +{ + int32_t h, e; +} eh_t; + +int ksw_extend2_origin(int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int is_left, // 是不是向左扩展 + int m, // 碱基种类 (5) + const int8_t *mat, // 每个位置的query和target的匹配得分 m*m + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数 + int w, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, + int zdrop, + int h0, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 +{ + eh_t *eh; // score array + 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; + uint8_t *qmem, *ref, *seq; + assert(h0 > 0); + // allocate memory + qp = malloc(qlen * m); + eh = calloc(qlen + 1, 8); + qmem = malloc(qlen + tlen); + seq = (uint8_t *)&qmem[0]; + ref = (uint8_t *)&qmem[qlen]; + if (is_left) + { + for (i = 0; i < qlen; ++i) + seq[i] = query[qlen - 1 - i]; + for (i = 0; i < tlen; ++i) + ref[i] = target[tlen - 1 - i]; + } + else + { + for (i = 0; i < qlen; ++i) + seq[i] = query[i]; + for (i = 0; i < tlen; ++i) + ref[i] = target[i]; + } + // generate the query profile + for (k = i = 0; k < m; ++k) + { + const int8_t *p = &mat[k * m]; + for (j = 0; j < qlen; ++j) + qp[i++] = p[seq[j]]; + } + // fill the first row + eh[0].h = h0; + eh[1].h = h0 > oe_ins ? h0 - oe_ins : 0; + for (j = 2; j <= qlen && eh[j - 1].h > e_ins; ++j) + eh[j].h = eh[j - 1].h - e_ins; + // adjust $w if it is too large + k = m * m; + for (i = 0, max = 0; i < k; ++i) // get the max score + max = max > mat[i] ? max : mat[i]; + max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.); + max_ins = max_ins > 1 ? max_ins : 1; + w = w < max_ins ? w : max_ins; + max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.); + max_del = max_del > 1 ? max_del : 1; + w = w < max_del ? w : max_del; // TODO: is this necessary? + // printf("%d\n", w); + // DP loop + max = h0, max_i = max_j = -1; + max_ie = -1, gscore = -1; + max_off = 0; + beg = 0, end = qlen; + + for (i = 0; LIKELY(i < tlen); ++i) + { + int t, f = 0, h1, m = 0, mj = -1; + int8_t *q = &qp[ref[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) end = qlen; 没用 + // compute the first column + if (beg == 0) + { + h1 = h0 - (o_del + e_del * (i + 1)); + if (h1 < 0) + h1 = 0; + } + else + h1 = 0; + for (j = beg; LIKELY(j < end); ++j) + { + // At the beginning of the loop: eh[j] = { H(i-1,j-1), E(i,j) }, f = F(i,j) and h1 = H(i,j-1) + // Similar to SSE2-SW, cells are computed in the following order: + // H(i,j) = max{H(i-1,j-1)+S(i,j), E(i,j), F(i,j)} + // E(i+1,j) = max{H(i,j)-gapo, E(i,j)} - gape + // F(i,j+1) = max{H(i,j)-gapo, F(i,j)} - gape + eh_t *p = &eh[j]; + int h, M = p->h, e = p->e; // get H(i-1,j-1) and E(i-1,j) + p->h = h1; // set H(i,j-1) for the next row + M = M ? M + q[j] : 0; // separating H and M to disallow a cigar like "100M3I3D20M" + h = M > e ? M : e; // e and f are guaranteed to be non-negative, so h>=0 even if M<0 + h = h > f ? h : f; + h1 = h; // save H(i,j) to h1 for the next column + mj = m > h ? mj : j; // record the position where max score is achieved + m = m > h ? m : h; // m is stored at eh[mj+1] + t = M - oe_del; + t = t > 0 ? t : 0; + e -= e_del; + e = e > t ? e : t; // computed E(i+1,j) + p->e = e; // save E(i+1,j) for the next row + t = M - oe_ins; + t = t > 0 ? t : 0; + f -= e_ins; + f = f > t ? f : t; // computed F(i,j+1) + } + eh[end].h = h1; + eh[end].e = 0; + if (j == qlen) + { + max_ie = gscore > h1 ? max_ie : i; + gscore = gscore > h1 ? gscore : h1; + } + if (m == 0) + break; + if (m > max) + { + max = m, max_i = i, max_j = mj; + max_off = max_off > abs(mj - i) ? max_off : abs(mj - i); + } + else if (zdrop > 0) + { + if (i - max_i > mj - max_j) + { + if (max - m - ((i - max_i) - (mj - max_j)) * e_del > zdrop) + break; + } + else + { + if (max - m - ((mj - max_j) - (i - max_i)) * e_ins > zdrop) + break; + } + } + // 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; + // beg = 0; end = qlen; // uncomment this line for debugging + } + + free(eh); + free(qp); + free(qmem); + if (_qle) + *_qle = max_j + 1; + if (_tle) + *_tle = max_i + 1; + if (_gtle) + *_gtle = max_ie + 1; + if (_gscore) + *_gscore = gscore; + if (_max_off) + *_max_off = max_off; + return max; +} diff --git a/ksw_avx2_u8.c b/ksw_avx2_u8.c new file mode 100644 index 0000000..aeb1554 --- /dev/null +++ b/ksw_avx2_u8.c @@ -0,0 +1,450 @@ +#include +#include +#include +#include +#include +#include +#include + +#ifdef __GNUC__ +#define LIKELY(x) __builtin_expect((x), 1) +#define UNLIKELY(x) __builtin_expect((x), 0) +#else +#define LIKELY(x) (x) +#define UNLIKELY(x) (x) +#endif + +#undef MAX +#undef MIN +#define MAX(x, y) ((x) > (y) ? (x) : (y)) +#define MIN(x, y) ((x) < (y) ? (x) : (y)) +#define SIMD_WIDTH 32 + +static const uint8_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { + {0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0}, + {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff}}; + +// static const uint8_t reverse_mask[SIMD_WIDTH] = {1, 0, 3, 2, 5, 4, 7, 6, 9, 8, 11, 10, 13, 12, 15, 14, 1, 0, 3, 2, 5, 4, 7, 6, 9, 8, 11, 10, 13, 12, 15, 14}; +static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 15, 14, 13, 12, 11, 10, 9, 8}; + +// const int permute_mask = _MM_SHUFFLE(0, 1, 2, 3); +#define permute_mask _MM_SHUFFLE(0, 1, 2, 3) +// 初始化变量 +#define SIMD_INIT \ + int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \ + __m256i zero_vec; \ + __m256i max_vec; \ + __m256i oe_del_vec; \ + __m256i oe_ins_vec; \ + __m256i e_del_vec; \ + __m256i e_ins_vec; \ + __m256i h_vec_mask[SIMD_WIDTH]; \ + __m256i reverse_mask_vec; \ + zero_vec = _mm256_setzero_si256(); \ + oe_del_vec = _mm256_set1_epi8(oe_del); \ + oe_ins_vec = _mm256_set1_epi8(oe_ins); \ + e_del_vec = _mm256_set1_epi8(e_del); \ + e_ins_vec = _mm256_set1_epi8(e_ins); \ + __m256i match_sc_vec = _mm256_set1_epi8(a); \ + __m256i mis_sc_vec = _mm256_set1_epi8(b); \ + __m256i amb_sc_vec = _mm256_set1_epi8(1); \ + __m256i amb_vec = _mm256_set1_epi8(4); \ + reverse_mask_vec = _mm256_loadu_si256((__m256i *)(reverse_mask)); \ + for (i = 0; i < SIMD_WIDTH; ++i) \ + h_vec_mask[i] = _mm256_loadu_si256((__m256i *)(&h_vec_int_mask[i])); + +/* + * e 表示当前ref的碱基被删除 + * f 表示当前seq的碱基插入 + * m 表示当前碱基匹配（可以相等，也可以不想等） + * h 表示最大值 + */ +// load向量化数据 +#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 - 1])); \ + __m256i ts_vec = _mm256_loadu_si256((__m256i *)(&ref[i])); + +// 比对ref和seq的序列，计算罚分 +#define SIMD_CMP_SEQ \ + ts_vec = _mm256_permute4x64_epi64(ts_vec, permute_mask); \ + ts_vec = _mm256_shuffle_epi8(ts_vec, reverse_mask_vec); \ + __m256i match_mask_vec = _mm256_cmpeq_epi8(qs_vec, ts_vec); \ + __m256i mis_score_vec = _mm256_andnot_si256(match_mask_vec, mis_sc_vec); \ + __m256i match_score_vec = _mm256_and_si256(match_sc_vec, match_mask_vec); \ + __m256i q_amb_mask_vec = _mm256_cmpeq_epi8(qs_vec, amb_vec); \ + __m256i t_amb_mask_vec = _mm256_cmpeq_epi8(ts_vec, amb_vec); \ + __m256i amb_mask_vec = _mm256_or_si256(q_amb_mask_vec, t_amb_mask_vec); \ + __m256i amb_score_vec = _mm256_and_si256(amb_mask_vec, amb_sc_vec); \ + mis_score_vec = _mm256_andnot_si256(amb_mask_vec, mis_score_vec); \ + mis_score_vec = _mm256_or_si256(amb_score_vec, mis_score_vec); \ + match_score_vec = _mm256_andnot_si256(amb_mask_vec, match_score_vec); + +// 向量化计算h, e, f, m +#define SIMD_COMPUTE \ + __m256i en_vec0 = _mm256_max_epu8(m1, oe_del_vec); \ + en_vec0 = _mm256_subs_epu8(en_vec0, oe_del_vec); \ + __m256i en_vec1 = _mm256_max_epu8(e1, e_del_vec); \ + en_vec1 = _mm256_subs_epu8(en_vec1, e_del_vec); \ + __m256i en_vec = _mm256_max_epu8(en_vec0, en_vec1); \ + __m256i fn_vec0 = _mm256_max_epu8(m1j1, oe_ins_vec); \ + fn_vec0 = _mm256_subs_epu8(fn_vec0, oe_ins_vec); \ + __m256i fn_vec1 = _mm256_max_epu8(f1j1, e_ins_vec); \ + fn_vec1 = _mm256_subs_epu8(fn_vec1, e_ins_vec); \ + __m256i fn_vec = _mm256_max_epu8(fn_vec0, fn_vec1); \ + __m256i mn_vec0 = _mm256_adds_epu8(h0j1, match_score_vec); \ + mn_vec0 = _mm256_max_epu8(mn_vec0, mis_score_vec); \ + mn_vec0 = _mm256_subs_epu8(mn_vec0, mis_score_vec); \ + __m256i mn_mask = _mm256_cmpeq_epi8(h0j1, zero_vec); \ + __m256i mn_vec = _mm256_andnot_si256(mn_mask, mn_vec0); \ + __m256i hn_vec0 = _mm256_max_epu8(en_vec, fn_vec); \ + __m256i hn_vec = _mm256_max_epu8(hn_vec0, mn_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[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 \ + uint8_t *maxVal = (uint8_t *)&max_vec; \ + max_vec = _mm256_max_epu8(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 1)); \ + max_vec = _mm256_max_epu8(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 2)); \ + max_vec = _mm256_max_epu8(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 3)); \ + max_vec = _mm256_max_epu8(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 4)); \ + max_vec = _mm256_max_epu8(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 5)); \ + max_vec = _mm256_max_epu8(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 6)); \ + max_vec = _mm256_max_epu8(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 7)); \ + max_vec = _mm256_max_epu8(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 8)); \ + max_vec = _mm256_max_epu8(max_vec, _mm256_permute2x128_si256(max_vec, max_vec, 0x01)); \ + m = maxVal[0]; \ + if (m > 0) \ + { \ + for (j = beg, i = iend; j <= end; j += SIMD_WIDTH, i -= SIMD_WIDTH) \ + { \ + __m256i h2_vec = _mm256_loadu_si256((__m256i *)(&hA2[j])); \ + __m256i vcmp = _mm256_cmpeq_epi8(h2_vec, max_vec); \ + uint32_t mask = _mm256_movemask_epi8(vcmp); \ + if (mask > 0) \ + { \ + int pos = SIMD_WIDTH - 1 - __builtin_clz(mask); \ + mj = j - 1 + pos; \ + mi = i - 1 - pos; \ + } \ + } \ + } + +// 每轮迭代后，交换数组 +#define SWAP_DATA_POINTER \ + uint8_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_avx2_u8(int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int is_left, // 是不是向左扩展 + int m, // 碱基种类 (5) + const int8_t *mat, // 每个位置的query和target的匹配得分 m*m + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数 + int a, // 碱基match时的分数 + int b, // 碱基mismatch时的惩罚分数（正数） + int w, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, + int zdrop, + int h0, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 +{ + uint8_t *mA, *hA, *eA, *fA, *mA1, *mA2, *hA0, *hA1, *eA1, *fA1, *hA2, *eA2, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M + 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, k, 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; + int val_mem_size = (col_size * 9 + 31) >> 5 << 5; // 32字节的整数倍 + int mem_size = seq_size + ref_size + val_mem_size; + + SIMD_INIT; // 初始化simd用的数据 + + assert(h0 > 0); + + // allocate memory + mem = malloc(mem_size); + qtmem = &mem[0]; + seq = (uint8_t *)&qtmem[0]; + ref = (uint8_t *)&qtmem[seq_size]; + if (is_left) + { + 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]; + } + else + { + for (i = 0; i < qlen; ++i) + seq[i] = query[i]; + for (i = 0; i < tlen; ++i) + ref[i + SIMD_WIDTH] = target[i]; + } + + vmem = &ref[ref_size]; + for (i = 0; i < val_mem_size; i += SIMD_WIDTH) + { + _mm256_storeu_si256((__m256i *)&vmem[i], zero_vec); + } + + hA = &vmem[0]; + mA = &vmem[col_size * 3]; + eA = &vmem[col_size * 5]; + fA = &vmem[col_size * 7]; + + hA0 = &hA[0]; + hA1 = &hA[col_size]; + hA2 = &hA1[col_size]; + mA1 = &mA[0]; + mA2 = &mA[col_size]; + eA1 = &eA[0]; + eA2 = &eA[col_size]; + fA1 = &fA[0]; + fA2 = &fA[col_size]; + + // adjust $w if it is too large + k = m * m; + // get the max score + for (i = 0, max = 0; i < k; ++i) + max = max > mat[i] ? max : mat[i]; + max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.); + max_ins = max_ins > 1 ? max_ins : 1; + w = w < max_ins ? w : max_ins; + max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.); + max_del = max_del > 1 ? max_del : 1; + w = w < max_del ? w : max_del; // TODO: is this necessary? + if (tlen < qlen) + w = MIN(tlen - 1, w); + + // DP loop + max = h0, max_i = max_j = -1; + max_ie = -1, gscore = -1; + ; + max_off = 0; + beg = 1; + end = qlen; + // init h0 + hA0[0] = h0; // 左上角 + + if (qlen == 0 || tlen == 0) + Dloop = 0; // 防止意外情况 + if (w >= qlen) + { + max_ie = 0; + gscore = 0; + } + + int m_last = 0; + int iend; + + for (D = 1; LIKELY(D < Dloop); ++D) + { + // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 + if (D > tlen) + { + span = MIN(Dloop - D, w); + beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1); + } + else + { + span = MIN(D - 1, w); + beg1 = MAX(1, ((D - w) / 2) + 1); + } + end1 = MIN(qlen, beg1 + span); + + if (beg < beg1) + beg = beg1; + if (end > end1) + end = end1; + if (beg > end) + break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug + + iend = D - (beg - 1); // ref开始计算的位置，倒序 + span = end - beg; + iStart = iend - span - 1; // 0开始的ref索引位置 + + // 每一轮需要记录的数据 + int m = 0, mj = -1, mi = -1; + max_vec = zero_vec; + + // 要处理边界 + // 左边界 处理f (insert) + if (iStart == 0) + { + hA1[end] = MAX(0, h0 - (o_ins + e_ins * end)); + } + // 上边界 + if (beg == 1) + { + hA1[0] = MAX(0, h0 - (o_del + e_del * iend)); + } + else + { + hA1[beg - 1] = 0; + eA1[beg - 1] = 0; + } + + for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH) + { + // 取数据 + SIMD_LOAD; + // 比对seq，计算罚分 + SIMD_CMP_SEQ; + // 计算 + SIMD_COMPUTE; + // 存储结果 + SIMD_STORE; + } + // 剩下的计算单元 + if (j <= end) + { + // 取数据 + SIMD_LOAD; + // 比对seq，计算罚分 + SIMD_CMP_SEQ; + // 计算 + SIMD_COMPUTE; + // 去除多余计算的部分 + SIMD_REMOVE_EXTRA; + // 存储结果 + SIMD_STORE; + } + + SIMD_FIND_MAX; + + // 注意最后跳出循环j的值 + j = end + 1; + + if (j == qlen + 1) + { + max_ie = gscore > hA2[qlen] ? max_ie : iStart; + gscore = gscore > hA2[qlen] ? gscore : hA2[qlen]; + } + if (m == 0 && m_last == 0) + break; // 一定要注意，斜对角遍历和按列遍历的不同点 + if (m > max) + { + max = m, max_i = mi, max_j = mj; + max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi); + } + else if (zdrop > 0) + { + if (mi - max_i > mj - max_j) + { + if (max - m - ((mi - max_i) - (mj - max_j)) * e_del > zdrop) + break; + } + else + { + if (max - m - ((mj - max_j) - (mi - max_i)) * e_ins > zdrop) + break; + } + } + + // 调整计算的边界 + 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; + + m_last = m; + // swap m, h, e, f + SWAP_DATA_POINTER; + } + + free(mem); + if (_qle) + *_qle = max_j + 1; + if (_tle) + *_tle = max_i + 1; + if (_gtle) + *_gtle = max_ie + 1; + if (_gscore) + *_gscore = gscore; + if (_max_off) + *_max_off = max_off; + return max; +} diff --git a/ksw_cuda.c b/ksw_cuda.c new file mode 100644 index 0000000..e69de29 diff --git a/ksw_normal.c b/ksw_normal.c new file mode 100644 index 0000000..aee34e6 --- /dev/null +++ b/ksw_normal.c @@ -0,0 +1,148 @@ +#include +#include +#include + +#ifdef __GNUC__ +#define LIKELY(x) __builtin_expect((x), 1) +#define UNLIKELY(x) __builtin_expect((x), 0) +#else +#define LIKELY(x) (x) +#define UNLIKELY(x) (x) +#endif + +typedef struct +{ + int32_t h, e; +} eh_t; + +int ksw_normal(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int o_del, int e_del, int o_ins, int e_ins, int w, int end_bonus, int zdrop, int h0, int *_qle, int *_tle, int *_gtle, int *_gscore, int *_max_off) +{ + eh_t *eh; // score array + 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); + // generate the query profile + for (k = i = 0; k < m; ++k) + { + const int8_t *p = &mat[k * m]; + for (j = 0; j < qlen; ++j) + qp[i++] = p[query[j]]; + } + // fill the first row + eh[0].h = h0; + eh[1].h = h0 > oe_ins ? h0 - oe_ins : 0; + for (j = 2; j <= qlen && eh[j - 1].h > e_ins; ++j) + eh[j].h = eh[j - 1].h - e_ins; + // adjust $w if it is too large + k = m * m; + for (i = 0, max = 0; i < k; ++i) // get the max score + max = max > mat[i] ? max : mat[i]; + max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.); + max_ins = max_ins > 1 ? max_ins : 1; + w = w < max_ins ? w : max_ins; + max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.); + max_del = max_del > 1 ? max_del : 1; + w = w < max_del ? w : max_del; // TODO: is this necessary? + // DP loop + max = h0, max_i = max_j = -1; + max_ie = -1, gscore = -1; + max_off = 0; + beg = 0, end = qlen; + for (i = 0; LIKELY(i < tlen); ++i) + { + int t, f = 0, h1, m = 0, mj = -1; + 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) + end = qlen; + // compute the first column + if (beg == 0) + { + h1 = h0 - (o_del + e_del * (i + 1)); + if (h1 < 0) + h1 = 0; + } + else + h1 = 0; + for (j = beg; LIKELY(j < end); ++j) + { + // At the beginning of the loop: eh[j] = { H(i-1,j-1), E(i,j) }, f = F(i,j) and h1 = H(i,j-1) + // Similar to SSE2-SW, cells are computed in the following order: + // H(i,j) = max{H(i-1,j-1)+S(i,j), E(i,j), F(i,j)} + // E(i+1,j) = max{H(i,j)-gapo, E(i,j)} - gape + // F(i,j+1) = max{H(i,j)-gapo, F(i,j)} - gape + eh_t *p = &eh[j]; + int h, M = p->h, e = p->e; // get H(i-1,j-1) and E(i-1,j) + p->h = h1; // set H(i,j-1) for the next row + M = M ? M + q[j] : 0; // separating H and M to disallow a cigar like "100M3I3D20M" + h = M > e ? M : e; // e and f are guaranteed to be non-negative, so h>=0 even if M<0 + h = h > f ? h : f; + h1 = h; // save H(i,j) to h1 for the next column + mj = m > h ? mj : j; // record the position where max score is achieved + m = m > h ? m : h; // m is stored at eh[mj+1] + t = M - oe_del; + t = t > 0 ? t : 0; + e -= e_del; + e = e > t ? e : t; // computed E(i+1,j) + p->e = e; // save E(i+1,j) for the next row + t = M - oe_ins; + t = t > 0 ? t : 0; + f -= e_ins; + f = f > t ? f : t; // computed F(i,j+1) + } + eh[end].h = h1; + eh[end].e = 0; + if (j == qlen) + { + max_ie = gscore > h1 ? max_ie : i; + gscore = gscore > h1 ? gscore : h1; + } + if (m == 0) + break; + if (m > max) + { + max = m, max_i = i, max_j = mj; + max_off = max_off > abs(mj - i) ? max_off : abs(mj - i); + } + else if (zdrop > 0) + { + if (i - max_i > mj - max_j) + { + if (max - m - ((i - max_i) - (mj - max_j)) * e_del > zdrop) + break; + } + else + { + if (max - m - ((mj - max_j) - (i - max_i)) * e_ins > zdrop) + break; + } + } + // 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; + // beg = 0; end = qlen; // uncomment this line for debugging + } + free(eh); + free(qp); + if (_qle) + *_qle = max_j + 1; + if (_tle) + *_tle = max_i + 1; + if (_gtle) + *_gtle = max_ie + 1; + if (_gscore) + *_gscore = gscore; + if (_max_off) + *_max_off = max_off; + return max; +} diff --git a/main.c b/main.c new file mode 100644 index 0000000..51ac5e2 --- /dev/null +++ b/main.c @@ -0,0 +1,262 @@ +#include +#include +#include +#include +#include +#include "sys/time.h" + +#define SW_NORMAL 0 +#define SW_AVX2 1 +#define SW_CUDA 2 +#define SW_ALL 3 + +#define BLOCK_BUF_SIZE 1048576 +#define READ_BUF_SIZE 2048 +#define SEQ_BUF_SIZE (BLOCK_BUF_SIZE + READ_BUF_SIZE) + +#ifdef SHOW_PERF +// 用来调试，计算感兴趣部分的运行时间 +// 获取当前毫秒数 +int64_t get_mseconds() +{ + struct timeval tv; + gettimeofday(&tv, NULL); + return (int64_t)1000 * (tv.tv_sec + ((1e-6) * tv.tv_usec)); +} + +int64_t time_sw_normal = 0, + time_sw_avx2 = 0, + time_sw_avx2_u8 = 0; + +#endif + +extern int ksw_normal(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int o_del, int e_del, int o_ins, int e_ins, int w, int end_bonus, int zdrop, int h0, int *_qle, int *_tle, int *_gtle, int *_gscore, int *_max_off); +extern int ksw_avx2(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int is_left, int m, const int8_t *mat, int o_del, int e_del, + int o_ins, int e_ins, int a, int b, int w, int end_bonus, int zdrop, int h0, int *_qle, int *_tle, int *_gtle, int *_gscore, int *_max_off); +extern int ksw_avx2_u8(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int is_left, int m, const int8_t *mat, int o_del, int e_del, + int o_ins, int e_ins, int a, int b, int w, int end_bonus, int zdrop, int h0, int *_qle, int *_tle, int *_gtle, int *_gscore, int *_max_off); + +/* + * 包含一个参数，用来区分调用那个sw算法 + * 参数为 normal/avx2/cuda + */ +// 程序执行入口 +int main(int argc, char *argv[]) +{ + /* + int sw_algo = SW_NORMAL; + + // 判断执行的sw的实现类型 + if (argc > 1) + { + if (strcmp(argv[1], "normal") == 0) + { + sw_algo = SW_NORMAL; + } + else if (strcmp(argv[1], "avx2") == 0) + { + sw_algo = SW_AVX2; + } + else if (strcmp(argv[1], "cuda") == 0) + { + sw_algo = SW_CUDA; + } + else + { + sw_algo = SW_ALL; + } + } */ + + // 初始化一些全局参数 + int8_t mat[25] = {1, -4, -4, -4, -1, + -4, 1, -4, -4, -1, + -4, -4, 1, -4, -1, + -4, -4, -4, 1, -1, + -1, -1, -1, -1, -1}; + int max_off[2]; + int qle, tle, gtle, gscore; + + // 读取测试数据 + char *query_arr = (char *)malloc(SEQ_BUF_SIZE); + char *target_arr = (char *)malloc(SEQ_BUF_SIZE); + int *info_buf = (int *)malloc(SEQ_BUF_SIZE); + int **info_arr = (int **)malloc(SEQ_BUF_SIZE); + FILE *query_f = 0, *target_f = 0, *info_f = 0; + // const char *qf_path = "/public/home/zzh/data/sw/q_s.fa"; + // const char *tf_path = "/public/home/zzh/data/sw/t_s.fa"; + // const char *if_path = "/public/home/zzh/data/sw/i_s.txt"; + const char *qf_path = "/public/home/zzh/data/sw/q_m.fa"; + const char *tf_path = "/public/home/zzh/data/sw/t_m.fa"; + const char *if_path = "/public/home/zzh/data/sw/i_m.txt"; + // const char *qf_path = "/public/home/zzh/data/sw/q_m.fa"; + // const char *tf_path = "/public/home/zzh/data/sw/t_m.fa"; + // const char *if_path = "/public/home/zzh/data/sw/i_m.txt"; + query_f = fopen(qf_path, "r"); + target_f = fopen(tf_path, "r"); + info_f = fopen(if_path, "r"); + + // 每次读取一定量的数据，然后执行，直到处理完所有数据 + int total_line_num = 0; // 目前处理的总的数据行数 + int block_line_num = 0; // 当前循环包含的数据行数 + int i, j; + // const int max_read = READ_BUF_SIZE; // 每次最多读取的字符 + char read_buf[READ_BUF_SIZE]; // 读文件缓存 + // int ret_code = 0; + + // 初始化info_arr数组 + i = 0; + j = 0; + while (1) + { + if (j > BLOCK_BUF_SIZE) + break; + info_arr[i] = &info_buf[j]; + i += 1; + j += 3; + } + + int score_normal = 0, score_avx2 = 0, score_avx2_u8 = 0; + + while (!feof(target_f)) + { + block_line_num = 0; + // target序列一般占用存储最多，先读取target，看一个buf能读多少行，query和info就按照这个行数来读 + int cur_read_size = 0; + while (!feof(target_f) && cur_read_size < BLOCK_BUF_SIZE) + { + if (fgets(read_buf, READ_BUF_SIZE, target_f) == NULL) + break; + const int line_size = strlen(read_buf); + assert(line_size < READ_BUF_SIZE); + ++block_line_num; + ++total_line_num; + strncpy(target_arr + cur_read_size, read_buf, line_size); + cur_read_size += line_size; + // fprintf(stderr, "%d %d \n", line_size, cur_read_size); + } + + // 读query + cur_read_size = 0; + for (i = 0; i < block_line_num; ++i) + { + if (fgets(read_buf, READ_BUF_SIZE, query_f) == NULL) + break; + const int line_size = strlen(read_buf); + assert(line_size < READ_BUF_SIZE); + strncpy(query_arr + cur_read_size, read_buf, line_size); + cur_read_size += line_size; + } + + // 读info + cur_read_size = 0; + for (i = 0; i < block_line_num; ++i) + { + if (fgets(read_buf, READ_BUF_SIZE, info_f) == NULL) + break; + const int line_size = strlen(read_buf); + assert(line_size < READ_BUF_SIZE); + sscanf(read_buf, "%d %d %d\n", &info_arr[i][0], &info_arr[i][1], &info_arr[i][2]); + cur_read_size += line_size; + // fprintf(stderr, "%-8d%-8d%-8d\n", info_arr[i][0], info_arr[i][1], info_arr[i][2]); + // fprintf(stderr, "%s\n", read_buf); + } + + // 性能测试 + + // 普通 sw + int cur_query_pos = 0; + int cur_target_pos = 0; + for (i = 0; i < block_line_num; ++i) + { +#ifdef SHOW_PERF + int64_t start_time = get_mseconds(); +#endif + score_normal += ksw_normal( + info_arr[i][0], + (uint8_t *)query_arr + cur_query_pos, + info_arr[i][1], + (uint8_t *)target_arr + cur_target_pos, + 5, mat, 6, 1, 6, 1, 100, 5, 100, + info_arr[i][2], + &qle, &tle, >le, &gscore, &max_off[0]); +#ifdef SHOW_PERF + time_sw_normal += get_mseconds() - start_time; +#endif + // 更新query和target位置信息 + cur_query_pos += info_arr[i][0]; + cur_target_pos += info_arr[i][1]; + // fprintf(stderr, "%d %d %d %d %d %d %d\n", score_normal, qle, tle, gtle, gscore, max_off[0], max_off[1]); + } + + // avx2 sw + cur_query_pos = 0; + cur_target_pos = 0; + for (i = 0; i < block_line_num; ++i) + { +#ifdef SHOW_PERF + int64_t start_time = get_mseconds(); +#endif + score_avx2 += ksw_avx2( + info_arr[i][0], + (uint8_t *)query_arr + cur_query_pos, + info_arr[i][1], + (uint8_t *)target_arr + cur_target_pos, + 0, 5, mat, 6, 1, 6, 1, + 1, 4, + 100, 5, 100, + info_arr[i][2], + &qle, &tle, >le, &gscore, &max_off[0]); +#ifdef SHOW_PERF + time_sw_avx2 += get_mseconds() - start_time; +#endif + // 更新query和target位置信息 + cur_query_pos += info_arr[i][0]; + cur_target_pos += info_arr[i][1]; + // fprintf(stderr, "%d %d %d %d %d %d %d\n", score_avx2, qle, tle, gtle, gscore, max_off[0], max_off[1]); + } + + // avx2 u8 sw + cur_query_pos = 0; + cur_target_pos = 0; + for (i = 0; i < block_line_num; ++i) + { +#ifdef SHOW_PERF + int64_t start_time = get_mseconds(); +#endif + score_avx2_u8 += ksw_avx2_u8( + info_arr[i][0], + (uint8_t *)query_arr + cur_query_pos, + info_arr[i][1], + (uint8_t *)target_arr + cur_target_pos, + 0, 5, mat, 6, 1, 6, 1, + 1, 4, + 100, 5, 100, + info_arr[i][2], + &qle, &tle, >le, &gscore, &max_off[0]); +#ifdef SHOW_PERF + time_sw_avx2_u8 += get_mseconds() - start_time; +#endif + // 更新query和target位置信息 + cur_query_pos += info_arr[i][0]; + cur_target_pos += info_arr[i][1]; + // fprintf(stderr, "%d %d %d %d %d %d %d\n", score_normal, qle, tle, gtle, gscore, max_off[0], max_off[1]); + } + + // fprintf(stderr, "%d %d \n", block_line_num, total_line_num); + } + + // fprintf(stderr, "%d \n", score_normal); + +#ifdef SHOW_PERF + fprintf(stderr, "time_sw_normal: %f s; score: %d\n", time_sw_normal / 1000.0, score_normal); + fprintf(stderr, "time_sw_avx2: %f s; score: %d\n", time_sw_avx2 / 1000.0, score_avx2); + fprintf(stderr, "time_sw_avx2_u8: %f s; score: %d\n", time_sw_avx2_u8 / 1000.0, score_avx2_u8); +#endif + + if (query_f != 0) + fclose(query_f); + if (target_f != 0) + fclose(target_f); + if (info_f != 0) + fclose(info_f); +} \ No newline at end of file