diff --git a/.vscode/settings.json b/.vscode/settings.json index 8854382..ec0fef7 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -6,6 +6,8 @@ "vector": "c", "__bit_reference": "c", "__split_buffer": "c", - "string": "c" + "string": "c", + "cstdint": "c", + "algorithm": "c" } } \ No newline at end of file diff --git a/Makefile b/Makefile index 8f00926..bfa68c9 100644 --- a/Makefile +++ b/Makefile @@ -2,13 +2,22 @@ CC= gcc #CFLAGS= -g -Wall -Wno-unused-function -mavx2 CFLAGS= -Wall -Wno-unused-function -O2 -mavx2 DFLAGS= -DSHOW_PERF -OBJS= ksw_normal.o ksw_avx2.o ksw_cuda.o ksw_avx2_u8.o bsw_avx2.o ksw_avx2_aligned.o thread_mem.o ksw_avx2_u8_aligned.o PROG= sw_perf -PROG2= sw_perf_discrete INCLUDES= LIBS= SUBDIRS= . +OBJS= ksw_ext_normal.o \ + ksw_ext_avx2.o \ + ksw_ext_avx2_u8.o \ + ksw_ext_cuda.o \ + ksw_ext_avx2_heuristics.o \ + ksw_ext_avx2_u8_heuristics.o \ + ksw_ext_avx2_aligned.o \ + ksw_ext_avx2_u8_aligned.o \ + thread_mem.o \ + utils.o + ifeq ($(shell uname -s),Linux) LIBS += -lrt endif @@ -18,16 +27,13 @@ endif .c.o: $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) $(CPPFLAGS) $< -o $@ -all:$(PROG) $(PROG2) +all:$(PROG) sw_perf:$(OBJS) main.o $(CC) $(CFLAGS) $(LDFLAGS) $(OBJS) main.o -o $@ -L. $(LIBS) -$(PROG2):$(OBJS) discrete_mem_main.o - $(CC) $(CFLAGS) $(LDFLAGS) $(OBJS) main.o -o $@ -L. $(LIBS) - clean: - rm -f *.o a.out $(PROG) $(PROG2) *~ *.a + rm -f *.o a.out $(PROG) *~ *.a depend: ( LC_ALL=C ; export LC_ALL; makedepend -Y -- $(CFLAGS) $(DFLAGS) $(CPPFLAGS) -- *.c ) diff --git a/discrete_mem_main.c b/discrete_mem_main.c deleted file mode 100644 index 9298ebf..0000000 --- a/discrete_mem_main.c +++ /dev/null @@ -1,243 +0,0 @@ -#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) - -// 将文件读取到离散的内存里，看一下对sw性能的影响 - -#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[]) -{ - - // 初始化一些全局参数 - 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); - target_arr[block_line_num] = (char *)malloc(line_size); - assert(line_size < READ_BUF_SIZE); - ++block_line_num; - ++total_line_num; - strncpy(target_arr[block_line_num], 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); - query_arr[i] = (char *)malloc(line_size); - strncpy(query_arr[i], 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[i], - info_arr[i][1], - (uint8_t *)target_arr[i], - 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[i], - info_arr[i][1], - (uint8_t *)target_arr[i], - 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[i], - info_arr[i][1], - (uint8_t *)target_arr[i], - 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 diff --git a/ksw_ext.h b/ksw_ext.h new file mode 100644 index 0000000..9f931d8 --- /dev/null +++ b/ksw_ext.h @@ -0,0 +1,164 @@ +/********************************************************************************************* + Description: Declarations of sw extend functions + + Copyright : All right reserved by NCIC.ICT + + Author : Zhang Zhonghai + Date : 2023/08/23 +***********************************************************************************************/ +#ifndef __KSW_EXT_H +#define __KSW_EXT_H +#include + +typedef struct _thread_mem_t thread_mem_t; + +// declaration of ksw functions + +int ksw_extend_normal(thread_mem_t *tmem, // 内存池 + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + 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 w, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int zdrop, // 如果比对过程中，太多mismatch，提前结束比对 + int h0, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off); // 取得最大得分时在query和reference上位置差的 最大值 + +int ksw_extend_avx2(thread_mem_t *tmem, // 内存池 + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off); // 取得最大得分时在query和reference上位置差的 最大值 + +int ksw_extend_avx2_u8(thread_mem_t *tmem, // 内存池 + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off); // 取得最大得分时在query和reference上位置差的 最大值 + +int ksw_extend_avx2_heuristics(thread_mem_t *tmem, // 内存池 + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off); // 取得最大得分时在query和reference上位置差的 最大值 + +int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, // 内存池 + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off); // 取得最大得分时在query和reference上位置差的 最大值 + +int ksw_extend_avx2_aligned(thread_mem_t *tmem, // 内存池 + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off); // 取得最大得分时在query和reference上位置差的 最大值 + +int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, // 内存池 + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off); // 取得最大得分时在query和reference上位置差的 最大值 + +#endif \ No newline at end of file diff --git a/ksw_avx2.c b/ksw_ext_avx2.c similarity index 54% rename from ksw_avx2.c rename to ksw_ext_avx2.c index 6c9482f..a938537 100644 --- a/ksw_avx2.c +++ b/ksw_ext_avx2.c @@ -5,6 +5,7 @@ #include #include #include +#include "thread_mem.h" #ifdef __GNUC__ #define LIKELY(x) __builtin_expect((x), 1) @@ -20,9 +21,6 @@ #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}, @@ -41,28 +39,28 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { {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 permute_mask _MM_SHUFFLE(0, 1, 2, 3) +#define permute_mask 27 // 初始化变量 -#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) \ +#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(base_match_score); \ + __m256i mis_sc_vec = _mm256_set1_epi16(-base_mis_score); \ + __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])); /* @@ -171,35 +169,33 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { 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上位置差的 最大值 +int ksw_extend_avx2(thread_mem_t *tmem, + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 { 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 i, iStart, 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; @@ -208,14 +204,15 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query SIMD_INIT; // 初始化simd用的数据 - assert(h0 > 0); + assert(init_score > 0); // allocate memory - mem = malloc(mem_size); + // mem = malloc(mem_size); + mem = thread_mem_request(tmem, mem_size); qtmem = (int16_t *)&mem[0]; seq = &qtmem[0]; ref = &qtmem[seq_size]; - if (is_left) + if (extend_left) { for (i = 0; i < qlen; ++i) seq[i] = query[qlen - 1 - i]; @@ -250,39 +247,37 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query fA1 = &fA[0]; fA2 = &fA[col_size]; - // adjust $w if it is too large - k = m * m; + // adjust $window_size if it is too large + // get the max score - for (i = 0, max = 0; i < k; ++i) - max = max > mat[i] ? max : mat[i]; + max = base_match_score; 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; + window_size = window_size < max_ins ? window_size : 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? + window_size = window_size < max_del ? window_size : max_del; // TODO: is this necessary? if (tlen < qlen) - w = MIN(tlen - 1, w); + window_size = MIN(tlen - 1, window_size); // DP loop - max = h0, max_i = max_j = -1; + max = init_score, max_i = max_j = -1; max_ie = -1, gscore = -1; ; max_off = 0; beg = 1; end = qlen; - // init h0 - hA0[0] = h0; // 左上角 + // init init_score + hA0[0] = init_score; // 左上角 if (qlen == 0 || tlen == 0) Dloop = 0; // 防止意外情况 - if (w >= qlen) + if (window_size >= qlen) { max_ie = 0; gscore = 0; } - int m_last = 0; int iend; for (D = 1; LIKELY(D < Dloop); ++D) @@ -290,13 +285,13 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 if (D > tlen) { - span = MIN(Dloop - D, w); - beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1); + span = MIN(Dloop - D, window_size); + beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1); } else { - span = MIN(D - 1, w); - beg1 = MAX(1, ((D - w) / 2) + 1); + span = MIN(D - 1, window_size); + beg1 = MAX(1, ((D - window_size) / 2) + 1); } end1 = MIN(qlen, beg1 + span); @@ -321,12 +316,12 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query // 左边界 处理f (insert) if (iStart == 0) { - hA1[end] = MAX(0, h0 - (o_ins + e_ins * end)); + hA1[end] = MAX(0, init_score - (o_ins + e_ins * end)); } // 上边界 if (beg == 1) { - hA1[0] = MAX(0, h0 - (o_del + e_del * iend)); + hA1[0] = MAX(0, init_score - (o_del + e_del * iend)); } else { @@ -370,26 +365,12 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query 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 (0) //(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) @@ -411,12 +392,12 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query // beg = 0; // end = qlen; // uncomment this line for debugging */ - m_last = m; // swap m, h, e, f SWAP_DATA_POINTER; } - free(mem); + // free(mem); + thread_mem_release(tmem, mem_size); if (_qle) *_qle = max_j + 1; if (_tle) @@ -428,182 +409,4 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query 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; -} +} \ No newline at end of file diff --git a/ksw_avx2_aligned.c b/ksw_ext_avx2_aligned.c similarity index 88% rename from ksw_avx2_aligned.c rename to ksw_ext_avx2_aligned.c index bbf9765..cd30a57 100644 --- a/ksw_avx2_aligned.c +++ b/ksw_ext_avx2_aligned.c @@ -198,28 +198,26 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { cur_match_arr = next_match_arr; \ next_match_arr = tmp; -// uint8_t mem_addr[102400]; - -int ksw_avx2_aligned(thread_mem_t *tmem, - int qlen, // query length 待匹配段碱基的query长度 - const uint8_t *query, // read碱基序列 - int tlen, // target length reference的长度 - const uint8_t *target, // reference序列 - int extend_left, // 是不是向左扩展 - int o_del, // deletion 错配开始的惩罚系数 - int e_del, // deletion extension的惩罚系数 - int o_ins, // insertion 错配开始的惩罚系数 - int e_ins, // insertion extension的惩罚系数SIMD_BTYES - int base_match_score, // 碱基match时的分数 - int base_mis_score, // 碱基mismatch时的惩罚分数（正数） - int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 - int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 - int init_score, // 该seed的初始得分（完全匹配query的碱基数） - int *_qle, // 匹配得到全局最大得分的碱基在query的位置 - int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 - int *_gtle, // query全部匹配上的target的长度 - int *_gscore, // query的端到端匹配得分 - int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 +int ksw_extend_avx2_aligned(thread_mem_t *tmem, + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + 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, @@ -333,8 +331,8 @@ int ksw_avx2_aligned(thread_mem_t *tmem, if (read_start_pos > read_end_pos) break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug - // read_start_pos = 1; - // read_end_pos = qlen; + 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索引位置 @@ -411,26 +409,26 @@ int ksw_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 = 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; // swap m, h, e, f SWAP_DATA_POINTER; } diff --git a/bsw_avx2.c b/ksw_ext_avx2_heuristics.c similarity index 71% rename from bsw_avx2.c rename to ksw_ext_avx2_heuristics.c index d7d846a..41ec171 100644 --- a/bsw_avx2.c +++ b/ksw_ext_avx2_heuristics.c @@ -5,6 +5,7 @@ #include #include #include +#include "thread_mem.h" #ifdef __GNUC__ #define LIKELY(x) __builtin_expect((x), 1) @@ -119,16 +120,6 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { hn_vec = _mm256_max_epi16(hn_vec, zero_vec); // 存储向量化结果 -// #define SIMD_STORE - -// __m256i cur_max_vec = _mm256_max_epu8(max_vec, hn_vec); \ -// __m256i vcmp = _mm256_cmpgt_epi8(cur_max_vec, max_vec); \ -// uint32_t mask = _mm256_movemask_epi8(vcmp); \ -// if (mask > 0) \ -// { \ -// simd_i = i; \ -// simd_j = j; \ -// } #define SIMD_STORE \ max_vec = _mm256_max_epu8(max_vec, hn_vec); \ _mm256_storeu_si256((__m256i *)&next_del_arr[j], en_vec); \ @@ -199,45 +190,45 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { cur_match_arr = next_match_arr; \ next_match_arr = tmp; -uint8_t mem[102400]; +// uint8_t mem[102400]; -int bsw_avx2(int qlen, // query length 待匹配段碱基的query长度 - const uint8_t *query, // read碱基序列 - int tlen, // target length reference的长度 - const uint8_t *target, // reference序列 - int extend_left, // 是不是向左扩展 - int o_del, // deletion 错配开始的惩罚系数 - int e_del, // deletion extension的惩罚系数 - int o_ins, // insertion 错配开始的惩罚系数 - int e_ins, // insertion extension的惩罚系数SIMD_BTYES - int base_match_score, // 碱基match时的分数 - int base_mis_score, // 碱基mismatch时的惩罚分数（正数） - int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 - int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 - int init_score, // 该seed的初始得分（完全匹配query的碱基数） - int *_qle, // 匹配得到全局最大得分的碱基在query的位置 - int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 - int *_gtle, // query全部匹配上的target的长度 - int *_gscore, // query的端到端匹配得分 - int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 +int ksw_extend_avx2_heuristics(thread_mem_t *tmem, + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 { - return 0; + // return 0; int16_t *mA, *eA, *hA, *fA, *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 *seq, *ref; - // uint8_t *mem; + 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 i, iStart, 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; int val_mem_size = (col_size * 9 * 2 + 31) >> 5 << 5; // 32字节的整数倍 int mem_size = (seq_size + ref_size) * 2 + val_mem_size; - int simd_i = -1, simd_j = -1; SIMD_INIT; // 初始化simd用的数据 @@ -245,6 +236,7 @@ int bsw_avx2(int qlen, // query length 待匹配段碱基的query // allocate memory // mem = malloc(mem_size); + mem = thread_mem_request(tmem, mem_size); qtmem = (int16_t *)&mem[0]; seq = &qtmem[0]; ref = &qtmem[seq_size]; @@ -313,7 +305,6 @@ int bsw_avx2(int qlen, // query length 待匹配段碱基的query gscore = 0; } - int m_last = 0; int iend; for (D = 1; LIKELY(D < Dloop); ++D) @@ -331,8 +322,8 @@ int bsw_avx2(int qlen, // query length 待匹配段碱基的query } end1 = MIN(qlen, beg1 + span); - // beg = 1; - // end = qlen; + beg = 1; + end = qlen; if (beg < beg1) beg = beg1; @@ -402,134 +393,38 @@ int bsw_avx2(int qlen, // query length 待匹配段碱基的query max_ie = gscore > next_max_arr[qlen] ? max_ie : iStart; gscore = gscore > next_max_arr[qlen] ? gscore : next_max_arr[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); } - // 调整计算的边界 - - // fprintf(stderr, "beg: %d, end: %d ", beg, end); - - /* for (j = beg; j <= end; j += SIMD_WIDTH) - { - __m256i h1 = _mm256_loadu_si256((__m256i *)(&cur_max_arr[j - 1])); - __m256i h2 = _mm256_loadu_si256((__m256i *)(&next_max_arr[j])); - __m256i orvec = _mm256_or_si256(h1, h2); - __m256i vcmp = _mm256_cmpgt_epi16(orvec, zero_vec); - uint32_t mask = _mm256_movemask_epi8(vcmp); - if (mask > 0) - { - // vcmp = _mm256_permute4x64_epi64(vcmp, permute_mask); - // vcmp = _mm256_shufflelo_epi16(vcmp, permute_mask); - // vcmp = _mm256_shufflehi_epi16(vcmp, permute_mask); - // mask = _mm256_movemask_epi8(vcmp); - //// int pos = SIMD_WIDTH - 1 - ((__builtin_clz(mask)) >> 1); - // int pos = ((__builtin_clz(mask)) >> 1); - // beg = j + pos; - int pos = __builtin_ctz(mask) >> 1; - beg = j + pos; - // if (beg > end) - // beg = end; - // beg = j + pos; - // beg = 0; - break; - } - } - */ - // beg = 0; - // fprintf(stderr, "new beg: %d\n", beg); - // int pos = 0; - // for (j = beg; j <= end; j += SIMD_WIDTH) + // for (j = beg; LIKELY(j <= end); ++j) //{ - // __m256i h1 = _mm256_loadu_si256((__m256i *)(&cur_max_arr[j - 1])); - // __m256i h2 = _mm256_loadu_si256((__m256i *)(&next_max_arr[j])); - // __m256i orvec = _mm256_or_si256(h1, h2); - // int *val = (int *)&orvec; - // for (i = 0; i < SIMD_WIDTH; ++i) - // if (val[i]) - // { - // pos = SIMD_WIDTH - 1 - i; - // break; - // } - //} + // int has_val = cur_max_arr[j - 1] | next_max_arr[j]; + // if (has_val) + // { + // break; + // } + // } // beg = j; - 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; - 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; - } - // else - // last_max_arr[j - 1] = 0; - } - end = j + 1 <= qlen ? j + 1 : qlen; - - /* for (j = end + 1; j >= beg; j -= SIMD_WIDTH) // 没有考虑beg附近，且长度小于SIMD_WIDTH的数据 - { - __m256i h1 = _mm256_loadu_si256((__m256i *)(&cur_max_arr[j - 1])); - __m256i h2 = _mm256_loadu_si256((__m256i *)(&next_max_arr[j])); - __m256i orvec = _mm256_or_si256(h1, h2); - __m256i vcmp = _mm256_cmpgt_epi16(orvec, zero_vec); - uint32_t mask = _mm256_movemask_epi8(vcmp); - if (mask > 0) - { - int pos = __builtin_clz(mask) >> 1; - const int new_end = j + SIMD_WIDTH - pos; - end = new_end <= qlen ? new_end : qlen; - break; - } - else - { - _mm256_storeu_si256((__m256i *)&last_max_arr[j - 1], zero_vec); - } - } - */ - // beg = 0; - // end = qlen; // uncomment this line for debugging - m_last = m; // swap m, h, e, f SWAP_DATA_POINTER; } - // __m256i origin_max_vec = max_vec; - // 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; - // max = maxVal[0]; - // if (max > 0) - // { - // __m256i vcmp = _mm256_cmpeq_epi16(origin_max_vec, max_vec); - // uint32_t mask = _mm256_movemask_epi8(vcmp); - // if (mask > 0) - // { - // int pos = SIMD_WIDTH - 1 - ((__builtin_clz(mask)) >> 1); - // int mj = simd_j - 1 + pos; - // int mi = simd_i - 1 - pos; - // max_i = mi, max_j = mj; - // max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi); - // } - // } - // free(mem); + thread_mem_release(tmem, mem_size); if (_qle) *_qle = max_j + 1; if (_tle) diff --git a/ksw_ext_avx2_u8.c b/ksw_ext_avx2_u8.c new file mode 100644 index 0000000..aac8e31 --- /dev/null +++ b/ksw_ext_avx2_u8.c @@ -0,0 +1,434 @@ +#include +#include +#include +#include +#include +#include +#include +#include "thread_mem.h" + +#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(base_match_score); \ + __m256i mis_sc_vec = _mm256_set1_epi8(base_mis_score); \ + __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_extend_avx2_u8(thread_mem_t *tmem, + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该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, 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(init_score > 0); + + // allocate memory + // mem = malloc(mem_size); + mem = thread_mem_request(tmem, mem_size); + qtmem = &mem[0]; + seq = (uint8_t *)&qtmem[0]; + ref = (uint8_t *)&qtmem[seq_size]; + if (extend_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 $window_size if it is too large + // get the max score + max = base_match_score; + max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.); + max_ins = max_ins > 1 ? max_ins : 1; + window_size = window_size < max_ins ? window_size : max_ins; + max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.); + max_del = max_del > 1 ? max_del : 1; + window_size = window_size < max_del ? window_size : max_del; // TODO: is this necessary? + if (tlen < qlen) + window_size = MIN(tlen - 1, window_size); + + // 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; // 左上角 + + if (qlen == 0 || tlen == 0) + Dloop = 0; // 防止意外情况 + if (window_size >= qlen) + { + max_ie = 0; + 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); + } + else + { + span = MIN(D - 1, window_size); + beg1 = MAX(1, ((D - window_size) / 2) + 1); + } + end1 = MIN(qlen, beg1 + span); + + if (beg < beg1) + beg = beg1; + if (end > end1) + end = end1; + if (beg > end) + break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug + + beg = 1; + end = qlen; + 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, 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; + } + + 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 > max) + { + max = m, max_i = mi, max_j = mj; + max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi); + } + + // 调整计算的边界 + // 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; + } + + // free(mem); + thread_mem_release(tmem, mem_size); + 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_aligned.c b/ksw_ext_avx2_u8_aligned.c similarity index 83% rename from ksw_avx2_u8_aligned.c rename to ksw_ext_avx2_u8_aligned.c index 72f5427..e4e140c 100644 --- a/ksw_avx2_u8_aligned.c +++ b/ksw_ext_avx2_u8_aligned.c @@ -69,17 +69,15 @@ static const uint8_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = { {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 permute_mask 27 // 初始化变量 #define SIMD_INIT \ int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \ __m256i zero_vec; \ - __m256i max_vec, last_max_vec = _mm256_set1_epi8(h0); \ + __m256i max_vec, last_max_vec = _mm256_set1_epi8(init_score); \ __m256i oe_del_vec; \ __m256i oe_ins_vec; \ __m256i e_del_vec; \ @@ -91,8 +89,8 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14, 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 match_sc_vec = _mm256_set1_epi8(base_match_score); \ + __m256i mis_sc_vec = _mm256_set1_epi8(base_mis_score); \ __m256i amb_sc_vec = _mm256_set1_epi8(1); \ __m256i amb_vec = _mm256_set1_epi8(4); \ reverse_mask_vec = _mm256_loadu_si256((__m256i *)(reverse_mask)); \ @@ -200,7 +198,7 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14, { \ if (read_seq[mj + 2] == ref_seq[mi + 1 + SIMD_WIDTH]) \ { \ - m += a; \ + m += base_match_score; \ } \ else \ { \ @@ -228,37 +226,31 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14, mA1 = mA2; \ mA2 = tmp; -int ksw_avx2_u8_aligned(thread_mem_t *tmem, - 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上位置差的 最大值 +int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem, + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该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 *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, k, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; + int i, iStart, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off; 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; uint8_t *mem_addr; int read_size = align_number(qlen * BASE_BYTES + MEM_ALIGN_BYTES); @@ -271,14 +263,14 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem, SIMD_INIT; // 初始化simd用的数据 - assert(h0 > 0); + assert(init_score > 0); mem_addr = thread_mem_request(tmem, request_mem_size); mem_addr = (void *)align_mem((uint64_t)mem_addr); ref_seq = (uint8_t *)&mem_addr[0]; read_seq = (uint8_t *)(mem_addr + ref_size + SIMD_BYTES - BASE_BYTES); - if (is_left) + if (extend_left) { for (i = 0; i < qlen; ++i) read_seq[i + 1] = query[qlen - 1 - i]; @@ -310,39 +302,36 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem, fA1 = &eA2[score_element_num]; fA2 = &fA1[score_element_num]; - // 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]; + // adjust $window_size if it is too large + + max = base_match_score; 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; + window_size = window_size < max_ins ? window_size : 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? + window_size = window_size < max_del ? window_size : max_del; // TODO: is this necessary? if (tlen < qlen) - w = MIN(tlen - 1, w); + window_size = MIN(tlen - 1, window_size); // DP loop - max = h0, max_i = max_j = -1; + max = init_score, max_i = max_j = -1; max_ie = -1, gscore = -1; ; max_off = 0; beg = 1; end = qlen; - // init h0 - hA0[0] = h0; // 左上角 + // init init_score + hA0[0] = init_score; // 左上角 if (qlen == 0 || tlen == 0) back_diagnal_num = 0; // 防止意外情况 - if (w >= qlen) + if (window_size >= qlen) { max_ie = 0; gscore = 0; } - int m_last = 0; int iend; for (D = 1; LIKELY(D < back_diagnal_num); ++D) @@ -350,13 +339,13 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem, // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 if (D > tlen) { - span = MIN(back_diagnal_num - D, w); - beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1); + span = MIN(back_diagnal_num - D, window_size); + beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1); } else { - span = MIN(D - 1, w); - beg1 = MAX(1, ((D - w) / 2) + 1); + span = MIN(D - 1, window_size); + beg1 = MAX(1, ((D - window_size) / 2) + 1); } end1 = MIN(qlen, beg1 + span); @@ -367,8 +356,8 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem, if (beg > end) break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug - // beg = 1; - // end = qlen; + beg = 1; + end = qlen; iend = D - (beg - 1); // ref开始计算的位置，倒序 span = end - beg; iStart = iend - span - 1; // 0开始的ref索引位置 @@ -381,12 +370,12 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem, // 左边界 处理f (insert) if (iStart == 0) { - hA1[end] = MAX(0, h0 - (o_ins + e_ins * end)); + hA1[end] = MAX(0, init_score - (o_ins + e_ins * end)); } // 上边界 if (beg == 1) { - hA1[0] = MAX(0, h0 - (o_del + e_del * iend)); + hA1[0] = MAX(0, init_score - (o_del + e_del * iend)); } else { @@ -438,26 +427,22 @@ int ksw_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; + // 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; - // beg = 0; - // end = qlen; - - m_last = m; // swap m, h, e, f SWAP_DATA_POINTER; } diff --git a/ksw_avx2_u8.c b/ksw_ext_avx2_u8_heuristics.c similarity index 82% rename from ksw_avx2_u8.c rename to ksw_ext_avx2_u8_heuristics.c index fb47ddd..9665a39 100644 --- a/ksw_avx2_u8.c +++ b/ksw_ext_avx2_u8_heuristics.c @@ -5,6 +5,7 @@ #include #include #include +#include "thread_mem.h" #ifdef __GNUC__ #define LIKELY(x) __builtin_expect((x), 1) @@ -64,7 +65,7 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14, #define SIMD_INIT \ int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \ __m256i zero_vec; \ - __m256i max_vec, last_max_vec = _mm256_set1_epi8(h0); \ + __m256i max_vec, last_max_vec = _mm256_set1_epi8(init_score); \ __m256i oe_del_vec; \ __m256i oe_ins_vec; \ __m256i e_del_vec; \ @@ -76,8 +77,8 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14, 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 match_sc_vec = _mm256_set1_epi8(base_match_score); \ + __m256i mis_sc_vec = _mm256_set1_epi8(base_mis_score); \ __m256i amb_sc_vec = _mm256_set1_epi8(1); \ __m256i amb_vec = _mm256_set1_epi8(4); \ reverse_mask_vec = _mm256_loadu_si256((__m256i *)(reverse_mask)); \ @@ -185,7 +186,7 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14, { \ if (seq[mj + 1] == ref[mi + 1 + SIMD_WIDTH]) \ { \ - m += a; \ + m += base_match_score; \ } \ else \ { \ @@ -213,45 +214,35 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14, mA1 = mA2; \ mA2 = tmp; -uint8_t mem1[102400]; +// uint8_t mem1[102400]; -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上位置差的 最大值 +int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem, + int qlen, // query length 待匹配段碱基的query长度 + const uint8_t *query, // read碱基序列 + int tlen, // target length reference的长度 + const uint8_t *target, // reference序列 + int extend_left, // 是不是向左扩展 + int o_del, // deletion 错配开始的惩罚系数 + int e_del, // deletion extension的惩罚系数 + int o_ins, // insertion 错配开始的惩罚系数 + int e_ins, // insertion extension的惩罚系数SIMD_BTYES + int base_match_score, // 碱基match时的分数 + int base_mis_score, // 碱基mismatch时的惩罚分数（正数） + int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离 + int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值 + int init_score, // 该seed的初始得分（完全匹配query的碱基数） + int *_qle, // 匹配得到全局最大得分的碱基在query的位置 + int *_tle, // 匹配得到全局最大得分的碱基在reference的位置 + int *_gtle, // query全部匹配上的target的长度 + int *_gscore, // query的端到端匹配得分 + int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值 { - // return h0; -#ifdef SHOW_PERF -// extern int64_t time_bsw_init; -// extern int64_t time_bsw_main_loop; -// extern int64_t time_bsw_find_max; -// extern int64_t time_bsw_adjust_bound; -// extern int64_t time_compare; -// int64_t start_time = get_mseconds(); -#endif + // return init_score; 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 = &mem1[0], *qtmem, *vmem; + 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 i, iStart, 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; @@ -260,14 +251,15 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que SIMD_INIT; // 初始化simd用的数据 - assert(h0 > 0); + assert(init_score > 0); // allocate memory // mem = malloc(mem_size); + mem = thread_mem_request(tmem, mem_size); qtmem = &mem[0]; seq = (uint8_t *)&qtmem[0]; ref = (uint8_t *)&qtmem[seq_size]; - if (is_left) + if (extend_left) { for (i = 0; i < qlen; ++i) seq[i] = query[qlen - 1 - i]; @@ -303,39 +295,36 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que fA1 = &fA[0]; fA2 = &fA[col_size]; - // adjust $w if it is too large - k = m * m; + // adjust $window_size if it is too large // get the max score - for (i = 0, max = 0; i < k; ++i) - max = max > mat[i] ? max : mat[i]; + max = base_match_score; 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; + window_size = window_size < max_ins ? window_size : 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? + window_size = window_size < max_del ? window_size : max_del; // TODO: is this necessary? if (tlen < qlen) - w = MIN(tlen - 1, w); + window_size = MIN(tlen - 1, window_size); // DP loop - max = h0, max_i = max_j = -1; + max = init_score, max_i = max_j = -1; max_ie = -1, gscore = -1; ; max_off = 0; beg = 1; end = qlen; - // init h0 - hA0[0] = h0; // 左上角 + // init init_score + hA0[0] = init_score; // 左上角 if (qlen == 0 || tlen == 0) Dloop = 0; // 防止意外情况 - if (w >= qlen) + if (window_size >= qlen) { max_ie = 0; gscore = 0; } - int m_last = 0; int iend; #ifdef SHOW_PERF // time_bsw_init += get_mseconds() - start_time; @@ -349,13 +338,13 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que // 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况 if (D > tlen) { - span = MIN(Dloop - D, w); - beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1); + span = MIN(Dloop - D, window_size); + beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1); } else { - span = MIN(D - 1, w); - beg1 = MAX(1, ((D - w) / 2) + 1); + span = MIN(D - 1, window_size); + beg1 = MAX(1, ((D - window_size) / 2) + 1); } end1 = MIN(qlen, beg1 + span); @@ -366,8 +355,8 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que if (beg > end) break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug - // beg = 1; - // end = qlen; + beg = 1; + end = qlen; iend = D - (beg - 1); // ref开始计算的位置，倒序 span = end - beg; iStart = iend - span - 1; // 0开始的ref索引位置 @@ -380,12 +369,12 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que // 左边界 处理f (insert) if (iStart == 0) { - hA1[end] = MAX(0, h0 - (o_ins + e_ins * end)); + hA1[end] = MAX(0, init_score - (o_ins + e_ins * end)); } // 上边界 if (beg == 1) { - hA1[0] = MAX(0, h0 - (o_del + e_del * iend)); + hA1[0] = MAX(0, init_score - (o_del + e_del * iend)); } else { @@ -418,19 +407,9 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que // 存储结果 SIMD_STORE; } -#ifdef SHOW_PERF -// time_bsw_main_loop += get_mseconds() - start_time; -#endif -#ifdef SHOW_PERF -// start_time = get_mseconds(); -#endif + SIMD_FIND_MAX; -#ifdef SHOW_PERF -// time_bsw_find_max += get_mseconds() - start_time; -#endif -#ifdef SHOW_PERF -// start_time = get_mseconds(); -#endif + // 注意最后跳出循环j的值 j = end + 1; @@ -439,50 +418,30 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que 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 (0) // (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; - hA2[end + 1] = 0; - 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; + // 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; - // beg = 0; - // end = qlen; - - m_last = m; // swap m, h, e, f SWAP_DATA_POINTER; #ifdef SHOW_PERF @@ -493,6 +452,7 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que } // free(mem); + thread_mem_release(tmem, mem_size); if (_qle) *_qle = max_j + 1; if (_tle) diff --git a/ksw_cuda.c b/ksw_ext_cuda.c similarity index 100% rename from ksw_cuda.c rename to ksw_ext_cuda.c diff --git a/ksw_normal.c b/ksw_ext_normal.c similarity index 92% rename from ksw_normal.c rename to ksw_ext_normal.c index 06a947c..5b612f6 100644 --- a/ksw_normal.c +++ b/ksw_ext_normal.c @@ -2,6 +2,7 @@ #include #include #include +#include "thread_mem.h" #ifdef __GNUC__ #define LIKELY(x) __builtin_expect((x), 1) @@ -16,15 +17,17 @@ 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) +int ksw_extend_normal(thread_mem_t *tmem, 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) { - return h0; + // return h0; 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); + // 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 for (k = i = 0; k < m; ++k) { @@ -140,8 +143,9 @@ int ksw_normal(int qlen, const uint8_t *query, int tlen, const uint8_t *target, // fprintf(stderr, "\n"); // fprintf(stderr, "%d\n", end); } - free(eh); - free(qp); + // free(eh); + // free(qp); + thread_mem_release(tmem, qlen * m + (qlen + 1) * 8); if (_qle) *_qle = max_j + 1; if (_tle) diff --git a/main.c b/main.c index 25164aa..d8dbd61 100644 --- a/main.c +++ b/main.c @@ -6,6 +6,8 @@ #include #include "sys/time.h" #include "thread_mem.h" +#include "ksw_ext.h" +#include "utils.h" #define SW_NORMAL 0 #define SW_AVX2 1 @@ -15,9 +17,12 @@ #define BLOCK_BUF_SIZE 1048576 #define READ_BUF_SIZE 2048 #define SEQ_BUF_SIZE (BLOCK_BUF_SIZE + READ_BUF_SIZE) +#define INIT_ALLOC_SIZE 4096 #define DIVIDE_BY (CLOCKS_PER_SEC * 1.0) +#define KERNEL_NUM 7 + #ifdef SHOW_PERF // 用来调试，计算感兴趣部分的运行时间 // 获取当前毫秒数 @@ -29,60 +34,89 @@ int64_t get_mseconds() return clock(); } -int64_t time_sw_normal = 0, - time_sw_avx2 = 0, - time_sw_avx2_u8 = 0, - time_sw_avx2_u8_aligned = 0, - time_bsw_avx2 = 0, - time_bsw_init = 0, - time_bsw_main_loop = 0, - time_bsw_find_max = 0, - time_bsw_adjust_bound = 0, - time_compare = 0; +int64_t time_sw[KERNEL_NUM] = {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); +#define _PERFORMANCE_TEST_NORMAL(kernel_num, func) \ + cur_query_pos = 0; \ + cur_target_pos = 0; \ + for (i = 0; i < block_line_num; ++i) \ + { \ + score[kernel_num] = func( \ + &tmem[kernel_num], \ + 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]); \ + score_total[kernel_num] += score[kernel_num]; \ + cur_query_pos += info_arr[i][0]; \ + cur_target_pos += info_arr[i][1]; \ + } -unsigned char nst_nt4_table[256] = { - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5 /*'-'*/, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4}; -char t_2bit2char[5] = {'A', 'C', 'G', 'T'}; -// 将碱基字符转成2位编码 -void convert_char_to_2bit(char *str) -{ - int i; - for (i = 0; i < strlen(str); ++i) - str[i] = nst_nt4_table[str[i]]; -} +#define _PERFORMANCE_TEST_AVX2(kernel_num, func) \ + cur_query_pos = 0; \ + cur_target_pos = 0; \ + for (i = 0; i < block_line_num; ++i) \ + { \ + score[kernel_num] = func( \ + &tmem[kernel_num], \ + info_arr[i][0], \ + (uint8_t *)query_arr + cur_query_pos, \ + info_arr[i][1], \ + (uint8_t *)target_arr + cur_target_pos, \ + 0, 6, 1, 6, 1, \ + 1, 4, \ + 100, 5, \ + info_arr[i][2], \ + &qle, &tle, >le, &gscore, &max_off[0]); \ + score_total[kernel_num] += score[kernel_num]; \ + cur_query_pos += info_arr[i][0]; \ + cur_target_pos += info_arr[i][1]; \ + } -// 读取测试数据 -int read_data() +#ifdef SHOW_PERF +#define PERFORMANCE_TEST_NORMAL(kernel_num, func) \ + start_time = get_mseconds(); \ + _PERFORMANCE_TEST_NORMAL(kernel_num, func); \ + time_sw[kernel_num] += get_mseconds() - start_time + +#define PERFORMANCE_TEST_AVX2(kernel_num, func) \ + start_time = get_mseconds(); \ + _PERFORMANCE_TEST_AVX2(kernel_num, func); \ + time_sw[kernel_num] += get_mseconds() - start_time +#else +#define PERFORMANCE_TEST_NORMAL(kernel_num, func) _PERFORMANCE_TEST_NORMAL(kernel_num, func) +#define PERFORMANCE_TEST_AVX2(kernel_num, func) _PERFORMANCE_TEST_AVX2(kernel_num, func) +#endif + +// 读取一行序列数据 +int read_seq_line(char *read_buf, FILE *f_ptr, char *out_arr) { - return 0; + if (fgets(read_buf, READ_BUF_SIZE, f_ptr) == NULL) + return 0; + int line_size = strlen(read_buf); + assert(line_size < READ_BUF_SIZE); + if (read_buf[line_size - 1] == '\n') + { + read_buf[line_size - 1] = '\0'; + line_size--; + } + convert_char_to_2bit(read_buf); + strncpy(out_arr, read_buf, line_size); + return line_size; } // 程序执行入口 int main(int argc, char *argv[]) { + const char *qf_path = argv[1]; + const char *tf_path = argv[2]; + const char *if_path = argv[3]; + // 初始化一些全局参数 int8_t mat[25] = {1, -4, -4, -4, -1, -4, 1, -4, -4, -1, @@ -91,52 +125,41 @@ int main(int argc, char *argv[]) -1, -1, -1, -1, -1}; int max_off[2]; int qle, tle, gtle, gscore; - thread_mem_t tmem, tmem_u8; - init_thread_mem(&tmem); - init_thread_mem(&tmem_u8); - // thread_mem_init_alloc(&tmem_u8, 10960); + thread_mem_t tmem[KERNEL_NUM]; + int i, j; + for (i = 0; i < KERNEL_NUM; ++i) + { + thread_mem_init_alloc(tmem + i, INIT_ALLOC_SIZE); + } + // 记录计算出的分数 + int score[KERNEL_NUM] = {0}; + int score_total[KERNEL_NUM] = {0}; // 读取测试数据 char *query_arr = (char *)malloc(SEQ_BUF_SIZE); char *target_arr = (char *)malloc(SEQ_BUF_SIZE); int *info_buf = (int *)malloc(SEQ_BUF_SIZE * sizeof(int)); int **info_arr = (int **)malloc(SEQ_BUF_SIZE * sizeof(int *)); + FILE *query_f = 0, *target_f = 0, *info_f = 0; - // const char *qf_path = "q.fa"; - // const char *tf_path = "t.fa"; - // const char *if_path = "i.txt"; - // const char *qf_path = "bug_q.fa"; - // const char *tf_path = "bug_t.fa"; - // const char *if_path = "bug_i.txt"; - // const char *qf_path = "/home/zzh/data/sw/q_s.fa"; - // const char *tf_path = "/home/zzh/data/sw/t_s.fa"; - // const char *if_path = "/home/zzh/data/sw/i_s.txt"; - // const char *qf_path = "/home/zzh/data/sw/q_m.fa"; - // const char *tf_path = "/home/zzh/data/sw/t_m.fa"; - // const char *if_path = "/home/zzh/data/sw/i_m.txt"; - const char *qf_path = "/home/zzh/data/sw/q_l.fa"; - const char *tf_path = "/home/zzh/data/sw/t_l.fa"; - const char *if_path = "/home/zzh/data/sw/i_l.txt"; - // const char *qf_path = "/home/zzh/data/sw/query.fa"; - // const char *tf_path = "/home/zzh/data/sw/target.fa"; - // const char *if_path = "/home/zzh/data/sw/info.txt"; + FILE *normal_out_f = 0, *avx2_out_f = 0, *avx2_u8_out_f = 0; + query_f = fopen(qf_path, "r"); target_f = fopen(tf_path, "r"); info_f = fopen(if_path, "r"); // 将每次比对的得分等信息写入文件，进行debug - FILE *normal_out_f = fopen("normal_out.txt", "w"); - FILE *avx2_out_f = fopen("avx2_out.txt", "w"); - FILE *avx2_u8_out_f = fopen("avx2_u8_out.txt", "w"); - FILE *bsw_avx2_out_f = fopen("bsw_avx2_out.txt", "w"); + // normal_out_f = fopen("normal_out.txt", "w"); + // avx2_out_f = fopen("avx2_out.txt", "w"); + // avx2_u8_out_f = fopen("avx2_u8_out.txt", "w"); // 每次读取一定量的数据，然后执行，直到处理完所有数据 int total_line_num = 0; // 目前处理的总的数据行数 int block_line_num = 0; // 当前循环包含的数据行数 - int i, j; - // const int max_read = READ_BUF_SIZE; // 每次最多读取的字符 + int cur_query_pos, cur_target_pos; + int64_t start_time; char read_buf[READ_BUF_SIZE]; // 读文件缓存 - // int ret_code = 0; + // 初始化info_arr数组 i = 0; j = 0; @@ -147,55 +170,30 @@ int main(int argc, char *argv[]) info_arr[i] = &info_buf[j]; i += 1; j += 3; - // fprintf(stderr, "%d\t%d\n", i, j); } - int score_normal = 0, score_avx2 = 0, score_avx2_u8 = 0, score_bsw_avx2 = 0; - int score_normal_total = 0, score_avx2_total = 0, score_avx2_u8_total = 0, score_bsw_avx2_total = 0; - int score_avx2_u8_aligned = 0, score_avx2_u8_aligned_total = 0; - while (!feof(target_f)) { - // fprintf(stderr, "debug\n"); - block_line_num = 0; + 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) + int line_size = read_seq_line(read_buf, target_f, target_arr + cur_read_size); + if (line_size == 0) break; - int line_size = strlen(read_buf); - assert(line_size < READ_BUF_SIZE); - if (read_buf[line_size - 1] == '\n') - { - read_buf[line_size - 1] = '\0'; - line_size--; - } - convert_char_to_2bit(read_buf); + cur_read_size += line_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); - // fprintf(stderr, "%d %d \n", block_line_num, total_line_num); } - // fprintf(stderr, "here\n"); // 读query cur_read_size = 0; for (i = 0; i < block_line_num; ++i) { - if (fgets(read_buf, READ_BUF_SIZE, query_f) == NULL) + int line_size = read_seq_line(read_buf, query_f, query_arr + cur_read_size); + if (line_size == 0) break; - int line_size = strlen(read_buf); - assert(line_size < READ_BUF_SIZE); - if (read_buf[line_size - 1] == '\n') - { - read_buf[line_size - 1] = '\0'; - line_size--; - } - convert_char_to_2bit(read_buf); - strncpy(query_arr + cur_read_size, read_buf, line_size); cur_read_size += line_size; } @@ -215,199 +213,42 @@ int main(int argc, char *argv[]) // 性能测试 - 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 - // 普通 sw - 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 - score_normal_total += score_normal; - // fprintf(normal_out_f, "%d %d\n", info_arr[i][2], score_normal); - // fprintf(stderr, "%d %d %d %d %d %d %d\n", info_arr[i][2], score_normal, qle, tle, gtle, gscore, max_off[0]); + // normal sw + PERFORMANCE_TEST_NORMAL(0, ksw_extend_normal); -// #ifdef SHOW_PERF -// start_time = get_mseconds(); -// #endif -// score_bsw_avx2 = ksw_avx2_aligned( -// &tmem, -// info_arr[i][0], -// (uint8_t *)query_arr + cur_query_pos, -// info_arr[i][1], -// (uint8_t *)target_arr + cur_target_pos, -// 0, 6, 1, 6, 1, -// 1, 4, -// 100, 5, -// info_arr[i][2], -// &qle, &tle, >le, &gscore, &max_off[0]); -// #ifdef SHOW_PERF -// time_bsw_avx2 += get_mseconds() - start_time; -// #endif -// score_bsw_avx2_total += score_bsw_avx2; -// fprintf(avx2_out_f, "%d %d\n", info_arr[i][2], score_avx2); -// fprintf(stderr, "%d %d %d %d %d %d %d\n", info_arr[i][2], score_bsw_avx2_total, qle, tle, gtle, gscore, max_off[0]); -/**/ -#ifdef SHOW_PERF - start_time = get_mseconds(); -#endif - score_avx2 = bsw_avx2( - info_arr[i][0], - (uint8_t *)query_arr + cur_query_pos, - info_arr[i][1], - (uint8_t *)target_arr + cur_target_pos, - 0, 6, 1, 6, 1, - 1, 4, - 100, 5, - info_arr[i][2], - &qle, &tle, >le, &gscore, &max_off[0]); - // 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 - score_avx2_total += score_avx2; - // fprintf(avx2_out_f, "%d %d\n", info_arr[i][2], score_avx2); - // fprintf(stderr, "%d %d %d %d %d %d %d\n", info_arr[i][2], score_avx2, qle, tle, gtle, gscore, max_off[0]); + // avx2 + PERFORMANCE_TEST_AVX2(1, ksw_extend_avx2); -#ifdef SHOW_PERF - 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 - score_avx2_u8_total += score_avx2_u8; - // fprintf(avx2_u8_out_f, "%d %d %d %d %d %d\n", score_avx2_u8, qle, tle, gtle, gscore, max_off[0]); + // avx2 u8 + PERFORMANCE_TEST_AVX2(2, ksw_extend_avx2_u8); -#ifdef SHOW_PERF - start_time = get_mseconds(); -#endif - score_avx2_u8_aligned = ksw_avx2_u8_aligned( - &tmem_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_aligned += get_mseconds() - start_time; -#endif - score_avx2_u8_aligned_total += score_avx2_u8_aligned; - // fprintf(avx2_u8_out_f, "%d %d %d %d %d %d\n", score_avx2_u8, qle, tle, gtle, gscore, max_off[0]); + // avx2 heuristics + PERFORMANCE_TEST_AVX2(3, ksw_extend_avx2_heuristics); - // 更新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 heuristics + PERFORMANCE_TEST_AVX2(4, ksw_extend_avx2_u8_heuristics); - // 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); + // avx2 mem aligned + PERFORMANCE_TEST_AVX2(5, ksw_extend_avx2_aligned); + + // avx2 u8 mem aligned + PERFORMANCE_TEST_AVX2(6, ksw_extend_avx2_u8_aligned); } - // fprintf(stderr, "%d \n", score_normal); - #ifdef SHOW_PERF - fprintf(stderr, "time_sw_normal: %f s; score: %d\n", time_sw_normal / DIVIDE_BY, score_normal_total); - fprintf(stderr, "time_bsw_avx2: %f s; score: %d\n", time_bsw_avx2 / DIVIDE_BY, score_bsw_avx2_total); - fprintf(stderr, "time_sw_avx2: %f s; score: %d\n", time_sw_avx2 / DIVIDE_BY, score_avx2_total); - fprintf(stderr, "time_sw_avx2_u8: %f s; score: %d\n", time_sw_avx2_u8 / DIVIDE_BY, score_avx2_u8_total); - fprintf(stderr, "time_sw_avx2_u8_aligned: %f s; score: %d\n", time_sw_avx2_u8_aligned / DIVIDE_BY, score_avx2_u8_aligned_total); - fprintf(stderr, "thread mem capacity: %d\t%d\n", tmem.capacity, tmem_u8.capacity); - fprintf(stderr, "time_bsw_init: %f s\n", time_bsw_init / DIVIDE_BY); - // fprintf(stderr, "time_bsw_main_loop: %f s\n", (time_bsw_main_loop) / DIVIDE_BY); - // fprintf(stderr, "time_bsw_find_max: %f s\n", (time_bsw_find_max) / DIVIDE_BY); - // fprintf(stderr, "time_bsw_adjust_bound: %f s\n", (time_bsw_adjust_bound) / DIVIDE_BY); - // fprintf(stderr, "time_bsw_main_loop: %f s\n", (time_bsw_main_loop - time_compare) / DIVIDE_BY); - // fprintf(stderr, "time_bsw_find_max: %f s\n", (time_bsw_find_max - time_compare) / DIVIDE_BY); - // fprintf(stderr, "time_bsw_adjust_bound: %f s\n", (time_bsw_adjust_bound - time_compare) / DIVIDE_BY); - // fprintf(stderr, "time_compare: %f s\n", time_compare / DIVIDE_BY); + char *kernel_names[7] = { + "normal", + "avx2", + "avx2_u8", + "avx2_heuristics", + "avx2_u8_heuristics", + "avx2_aligned", + "avx2_u8_aligned"}; + + for (i = 0; i < KERNEL_NUM; ++i) + { + fprintf(stderr, "[%18s] time: %9.6f s; score: %d\n", kernel_names[i], time_sw[i] / DIVIDE_BY, score_total[i]); + } #endif if (query_f != 0) @@ -422,6 +263,4 @@ int main(int argc, char *argv[]) fclose(avx2_u8_out_f); if (normal_out_f != 0) fclose(normal_out_f); - if (bsw_avx2_out_f != 0) - fclose(bsw_avx2_out_f); } diff --git a/run_all.sh b/run_all.sh new file mode 100755 index 0000000..312b1ae --- /dev/null +++ b/run_all.sh @@ -0,0 +1,2 @@ +#!/bin/bash +/home/zzh/work/sw_perf/sw_perf /home/zzh/data/sw/query.fa /home/zzh/data/sw/target.fa /home/zzh/data/sw/info.txt \ No newline at end of file diff --git a/run_l.sh b/run_l.sh new file mode 100755 index 0000000..02285a8 --- /dev/null +++ b/run_l.sh @@ -0,0 +1,2 @@ +#!/bin/bash +/home/zzh/work/sw_perf/sw_perf /home/zzh/data/sw/q_l.fa /home/zzh/data/sw/t_l.fa /home/zzh/data/sw/i_l.txt \ No newline at end of file diff --git a/run_m.sh b/run_m.sh new file mode 100755 index 0000000..b936788 --- /dev/null +++ b/run_m.sh @@ -0,0 +1,2 @@ +#!/bin/bash +/home/zzh/work/sw_perf/sw_perf /home/zzh/data/sw/q_m.fa /home/zzh/data/sw/t_m.fa /home/zzh/data/sw/i_m.txt \ No newline at end of file diff --git a/run_s.sh b/run_s.sh new file mode 100755 index 0000000..d31fc52 --- /dev/null +++ b/run_s.sh @@ -0,0 +1,2 @@ +#!/bin/bash +/home/zzh/work/sw_perf/sw_perf /home/zzh/data/sw/q_s.fa /home/zzh/data/sw/t_s.fa /home/zzh/data/sw/i_s.txt \ No newline at end of file diff --git a/thread_mem.c b/thread_mem.c index 1491f7b..3a5a930 100644 --- a/thread_mem.c +++ b/thread_mem.c @@ -9,6 +9,7 @@ #include "thread_mem.h" #include +#include // 创建 thread_mem_t *create_thread_mem() @@ -64,6 +65,13 @@ void *thread_mem_request(thread_mem_t *tmem, size_t byte_cnt) return ret_mem; } +void *thread_mem_request_and_clean(thread_mem_t *tmem, size_t byte_cnt) +{ + void *mem = thread_mem_request(tmem, byte_cnt); + memset(mem, 0, byte_cnt); + return mem; +} + // 将不用的内存归还给thread mem void thread_mem_release(thread_mem_t *tmem, size_t byte_cnt) { diff --git a/thread_mem.h b/thread_mem.h index 3035cea..4f441f1 100644 --- a/thread_mem.h +++ b/thread_mem.h @@ -16,7 +16,7 @@ #define MEM_ALIGN_BYTE 8 #define MEM_MOVE_BIT 3 -typedef struct +typedef struct _thread_mem_t { size_t occupied; // 已经占用的容量（字节数） 对齐的 size_t capacity; // 总容量（字节数） @@ -29,9 +29,13 @@ thread_mem_t *create_thread_mem(); void init_thread_mem(thread_mem_t *tmem); // 初始化并开辟一定量的内存 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); +// 请求内存并初始化为零 +void *thread_mem_request_and_clean(thread_mem_t *tmem, size_t byte_cnt); + // 将不用的内存归还给thread mem void thread_mem_release(thread_mem_t *tmem, size_t byte_cnt); diff --git a/utils.c b/utils.c new file mode 100644 index 0000000..33f0865 --- /dev/null +++ b/utils.c @@ -0,0 +1,39 @@ +/********************************************************************************************* + Description: Some useful functions + + Copyright : All right reserved by NCIC.ICT + + Author : Zhang Zhonghai + Date : 2023/08/25 +***********************************************************************************************/ +#include "utils.h" +#include +#include + +unsigned char nst_nt4_table[256] = { + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5 /*'-'*/, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4}; + +char t_2bit2char[5] = {'A', 'C', 'G', 'T', 'N'}; + +// 将碱基字符转成2位编码 +void convert_char_to_2bit(char *str) +{ + int i; + for (i = 0; i < strlen(str); ++i) + str[i] = nst_nt4_table[(uint8_t)str[i]]; +} \ No newline at end of file diff --git a/utils.h b/utils.h new file mode 100644 index 0000000..79ca7f7 --- /dev/null +++ b/utils.h @@ -0,0 +1,15 @@ +/********************************************************************************************* + Description: Some useful functions + + Copyright : All right reserved by NCIC.ICT + + Author : Zhang Zhonghai + Date : 2023/08/25 +***********************************************************************************************/ +#ifndef __UTILS_H +#define __UTILS_H + +// 将碱基字符转成2位编码 +void convert_char_to_2bit(char *str); + +#endif \ No newline at end of file