代码重构，总共7个kernel
This commit is contained in:
parent
b95e622e7a
commit
b53569db63
|
|
@ -6,6 +6,8 @@
|
||||||
"vector": "c",
|
"vector": "c",
|
||||||
"__bit_reference": "c",
|
"__bit_reference": "c",
|
||||||
"__split_buffer": "c",
|
"__split_buffer": "c",
|
||||||
"string": "c"
|
"string": "c",
|
||||||
|
"cstdint": "c",
|
||||||
|
"algorithm": "c"
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
20
Makefile
20
Makefile
|
|
@ -2,12 +2,21 @@ CC= gcc
|
||||||
#CFLAGS= -g -Wall -Wno-unused-function -mavx2
|
#CFLAGS= -g -Wall -Wno-unused-function -mavx2
|
||||||
CFLAGS= -Wall -Wno-unused-function -O2 -mavx2
|
CFLAGS= -Wall -Wno-unused-function -O2 -mavx2
|
||||||
DFLAGS= -DSHOW_PERF
|
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
|
PROG= sw_perf
|
||||||
PROG2= sw_perf_discrete
|
|
||||||
INCLUDES=
|
INCLUDES=
|
||||||
LIBS=
|
LIBS=
|
||||||
SUBDIRS= .
|
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)
|
ifeq ($(shell uname -s),Linux)
|
||||||
LIBS += -lrt
|
LIBS += -lrt
|
||||||
|
|
@ -18,16 +27,13 @@ endif
|
||||||
.c.o:
|
.c.o:
|
||||||
$(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) $(CPPFLAGS) $< -o $@
|
$(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) $(CPPFLAGS) $< -o $@
|
||||||
|
|
||||||
all:$(PROG) $(PROG2)
|
all:$(PROG)
|
||||||
|
|
||||||
sw_perf:$(OBJS) main.o
|
sw_perf:$(OBJS) main.o
|
||||||
$(CC) $(CFLAGS) $(LDFLAGS) $(OBJS) main.o -o $@ -L. $(LIBS)
|
$(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:
|
clean:
|
||||||
rm -f *.o a.out $(PROG) $(PROG2) *~ *.a
|
rm -f *.o a.out $(PROG) *~ *.a
|
||||||
|
|
||||||
depend:
|
depend:
|
||||||
( LC_ALL=C ; export LC_ALL; makedepend -Y -- $(CFLAGS) $(DFLAGS) $(CPPFLAGS) -- *.c )
|
( LC_ALL=C ; export LC_ALL; makedepend -Y -- $(CFLAGS) $(DFLAGS) $(CPPFLAGS) -- *.c )
|
||||||
|
|
|
||||||
|
|
@ -1,243 +0,0 @@
|
||||||
#include <stdlib.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <stdint.h>
|
|
||||||
#include <assert.h>
|
|
||||||
#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);
|
|
||||||
}
|
|
||||||
|
|
@ -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 <stdint.h>
|
||||||
|
|
||||||
|
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
|
||||||
|
|
@ -5,6 +5,7 @@
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include <immintrin.h>
|
#include <immintrin.h>
|
||||||
#include <emmintrin.h>
|
#include <emmintrin.h>
|
||||||
|
#include "thread_mem.h"
|
||||||
|
|
||||||
#ifdef __GNUC__
|
#ifdef __GNUC__
|
||||||
#define LIKELY(x) __builtin_expect((x), 1)
|
#define LIKELY(x) __builtin_expect((x), 1)
|
||||||
|
|
@ -20,9 +21,6 @@
|
||||||
#define MIN(x, y) ((x) < (y) ? (x) : (y))
|
#define MIN(x, y) ((x) < (y) ? (x) : (y))
|
||||||
#define SIMD_WIDTH 16
|
#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] = {
|
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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||||
{0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
{0xffff, 0xffff, 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, 0},
|
||||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}};
|
{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 \
|
#define SIMD_INIT \
|
||||||
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
||||||
__m256i zero_vec; \
|
__m256i zero_vec; \
|
||||||
__m256i max_vec; \
|
__m256i max_vec; \
|
||||||
__m256i oe_del_vec; \
|
__m256i oe_del_vec; \
|
||||||
__m256i oe_ins_vec; \
|
__m256i oe_ins_vec; \
|
||||||
__m256i e_del_vec; \
|
__m256i e_del_vec; \
|
||||||
__m256i e_ins_vec; \
|
__m256i e_ins_vec; \
|
||||||
__m256i h_vec_mask[SIMD_WIDTH]; \
|
__m256i h_vec_mask[SIMD_WIDTH]; \
|
||||||
zero_vec = _mm256_setzero_si256(); \
|
zero_vec = _mm256_setzero_si256(); \
|
||||||
oe_del_vec = _mm256_set1_epi16(-oe_del); \
|
oe_del_vec = _mm256_set1_epi16(-oe_del); \
|
||||||
oe_ins_vec = _mm256_set1_epi16(-oe_ins); \
|
oe_ins_vec = _mm256_set1_epi16(-oe_ins); \
|
||||||
e_del_vec = _mm256_set1_epi16(-e_del); \
|
e_del_vec = _mm256_set1_epi16(-e_del); \
|
||||||
e_ins_vec = _mm256_set1_epi16(-e_ins); \
|
e_ins_vec = _mm256_set1_epi16(-e_ins); \
|
||||||
__m256i match_sc_vec = _mm256_set1_epi16(a); \
|
__m256i match_sc_vec = _mm256_set1_epi16(base_match_score); \
|
||||||
__m256i mis_sc_vec = _mm256_set1_epi16(-b); \
|
__m256i mis_sc_vec = _mm256_set1_epi16(-base_mis_score); \
|
||||||
__m256i amb_sc_vec = _mm256_set1_epi16(-1); \
|
__m256i amb_sc_vec = _mm256_set1_epi16(-1); \
|
||||||
__m256i amb_vec = _mm256_set1_epi16(4); \
|
__m256i amb_vec = _mm256_set1_epi16(4); \
|
||||||
for (i = 0; i < SIMD_WIDTH; ++i) \
|
for (i = 0; i < SIMD_WIDTH; ++i) \
|
||||||
h_vec_mask[i] = _mm256_loadu_si256((__m256i *)(&h_vec_int_mask[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; \
|
mA1 = mA2; \
|
||||||
mA2 = tmp;
|
mA2 = tmp;
|
||||||
|
|
||||||
int ksw_avx2(int qlen, // query length 待匹配段碱基的query长度
|
int ksw_extend_avx2(thread_mem_t *tmem,
|
||||||
const uint8_t *query, // read碱基序列
|
int qlen, // query length 待匹配段碱基的query长度
|
||||||
int tlen, // target length reference的长度
|
const uint8_t *query, // read碱基序列
|
||||||
const uint8_t *target, // reference序列
|
int tlen, // target length reference的长度
|
||||||
int is_left, // 是不是向左扩展
|
const uint8_t *target, // reference序列
|
||||||
int m, // 碱基种类 (5)
|
int extend_left, // 是不是向左扩展
|
||||||
const int8_t *mat, // 每个位置的query和target的匹配得分 m*m
|
int o_del, // deletion 错配开始的惩罚系数
|
||||||
int o_del, // deletion 错配开始的惩罚系数
|
int e_del, // deletion extension的惩罚系数
|
||||||
int e_del, // deletion extension的惩罚系数
|
int o_ins, // insertion 错配开始的惩罚系数
|
||||||
int o_ins, // insertion 错配开始的惩罚系数
|
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||||
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
int base_match_score, // 碱基match时的分数
|
||||||
int a, // 碱基match时的分数
|
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||||
int b, // 碱基mismatch时的惩罚分数（正数）
|
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||||
int w, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||||
int end_bonus,
|
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||||
int zdrop,
|
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||||
int h0, // 该seed的初始得分（完全匹配query的碱基数）
|
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
int *_gtle, // query全部匹配上的target的长度
|
||||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
int *_gscore, // query的端到端匹配得分
|
||||||
int *_gtle, // query全部匹配上的target的长度
|
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
||||||
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 *mA, *hA, *eA, *fA, *mA1, *mA2, *hA0, *hA1, *eA1, *fA1, *hA2, *eA2, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M
|
||||||
int16_t *seq, *ref;
|
int16_t *seq, *ref;
|
||||||
uint8_t *mem;
|
uint8_t *mem;
|
||||||
int16_t *qtmem, *vmem;
|
int16_t *qtmem, *vmem;
|
||||||
int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH;
|
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 Dloop = tlen + qlen; // 循环跳出条件
|
||||||
int span, beg1, end1; // 边界条件计算
|
int span, beg1, end1; // 边界条件计算
|
||||||
int col_size = qlen + 2 + SIMD_WIDTH;
|
int col_size = qlen + 2 + SIMD_WIDTH;
|
||||||
|
|
@ -208,14 +204,15 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
|
|
||||||
SIMD_INIT; // 初始化simd用的数据
|
SIMD_INIT; // 初始化simd用的数据
|
||||||
|
|
||||||
assert(h0 > 0);
|
assert(init_score > 0);
|
||||||
|
|
||||||
// allocate memory
|
// allocate memory
|
||||||
mem = malloc(mem_size);
|
// mem = malloc(mem_size);
|
||||||
|
mem = thread_mem_request(tmem, mem_size);
|
||||||
qtmem = (int16_t *)&mem[0];
|
qtmem = (int16_t *)&mem[0];
|
||||||
seq = &qtmem[0];
|
seq = &qtmem[0];
|
||||||
ref = &qtmem[seq_size];
|
ref = &qtmem[seq_size];
|
||||||
if (is_left)
|
if (extend_left)
|
||||||
{
|
{
|
||||||
for (i = 0; i < qlen; ++i)
|
for (i = 0; i < qlen; ++i)
|
||||||
seq[i] = query[qlen - 1 - i];
|
seq[i] = query[qlen - 1 - i];
|
||||||
|
|
@ -250,39 +247,37 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
fA1 = &fA[0];
|
fA1 = &fA[0];
|
||||||
fA2 = &fA[col_size];
|
fA2 = &fA[col_size];
|
||||||
|
|
||||||
// adjust $w if it is too large
|
// adjust $window_size if it is too large
|
||||||
k = m * m;
|
|
||||||
// get the max score
|
// get the max score
|
||||||
for (i = 0, max = 0; i < k; ++i)
|
max = base_match_score;
|
||||||
max = max > mat[i] ? max : mat[i];
|
|
||||||
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.);
|
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.);
|
||||||
max_ins = max_ins > 1 ? max_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 = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.);
|
||||||
max_del = max_del > 1 ? max_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)
|
if (tlen < qlen)
|
||||||
w = MIN(tlen - 1, w);
|
window_size = MIN(tlen - 1, window_size);
|
||||||
|
|
||||||
// DP loop
|
// DP loop
|
||||||
max = h0, max_i = max_j = -1;
|
max = init_score, max_i = max_j = -1;
|
||||||
max_ie = -1, gscore = -1;
|
max_ie = -1, gscore = -1;
|
||||||
;
|
;
|
||||||
max_off = 0;
|
max_off = 0;
|
||||||
beg = 1;
|
beg = 1;
|
||||||
end = qlen;
|
end = qlen;
|
||||||
// init h0
|
// init init_score
|
||||||
hA0[0] = h0; // 左上角
|
hA0[0] = init_score; // 左上角
|
||||||
|
|
||||||
if (qlen == 0 || tlen == 0)
|
if (qlen == 0 || tlen == 0)
|
||||||
Dloop = 0; // 防止意外情况
|
Dloop = 0; // 防止意外情况
|
||||||
if (w >= qlen)
|
if (window_size >= qlen)
|
||||||
{
|
{
|
||||||
max_ie = 0;
|
max_ie = 0;
|
||||||
gscore = 0;
|
gscore = 0;
|
||||||
}
|
}
|
||||||
|
|
||||||
int m_last = 0;
|
|
||||||
int iend;
|
int iend;
|
||||||
|
|
||||||
for (D = 1; LIKELY(D < Dloop); ++D)
|
for (D = 1; LIKELY(D < Dloop); ++D)
|
||||||
|
|
@ -290,13 +285,13 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
||||||
if (D > tlen)
|
if (D > tlen)
|
||||||
{
|
{
|
||||||
span = MIN(Dloop - D, w);
|
span = MIN(Dloop - D, window_size);
|
||||||
beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1);
|
beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1);
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
span = MIN(D - 1, w);
|
span = MIN(D - 1, window_size);
|
||||||
beg1 = MAX(1, ((D - w) / 2) + 1);
|
beg1 = MAX(1, ((D - window_size) / 2) + 1);
|
||||||
}
|
}
|
||||||
end1 = MIN(qlen, beg1 + span);
|
end1 = MIN(qlen, beg1 + span);
|
||||||
|
|
||||||
|
|
@ -321,12 +316,12 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
// 左边界 处理f (insert)
|
// 左边界 处理f (insert)
|
||||||
if (iStart == 0)
|
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)
|
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
|
else
|
||||||
{
|
{
|
||||||
|
|
@ -370,26 +365,12 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
max_ie = gscore > hA2[qlen] ? max_ie : iStart;
|
max_ie = gscore > hA2[qlen] ? max_ie : iStart;
|
||||||
gscore = gscore > hA2[qlen] ? gscore : hA2[qlen];
|
gscore = gscore > hA2[qlen] ? gscore : hA2[qlen];
|
||||||
}
|
}
|
||||||
// if (m == 0 && m_last == 0)
|
|
||||||
// break; // 一定要注意，斜对角遍历和按列遍历的不同点
|
|
||||||
if (m > max)
|
if (m > max)
|
||||||
{
|
{
|
||||||
max = m, max_i = mi, max_j = mj;
|
max = m, max_i = mi, max_j = mj;
|
||||||
max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi);
|
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)
|
/*for (j = beg; LIKELY(j <= end); ++j)
|
||||||
|
|
@ -411,190 +392,12 @@ int ksw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
// beg = 0;
|
// beg = 0;
|
||||||
// end = qlen; // uncomment this line for debugging
|
// end = qlen; // uncomment this line for debugging
|
||||||
*/
|
*/
|
||||||
m_last = m;
|
|
||||||
// swap m, h, e, f
|
// swap m, h, e, f
|
||||||
SWAP_DATA_POINTER;
|
SWAP_DATA_POINTER;
|
||||||
}
|
}
|
||||||
|
|
||||||
free(mem);
|
// free(mem);
|
||||||
if (_qle)
|
thread_mem_release(tmem, mem_size);
|
||||||
*_qle = max_j + 1;
|
|
||||||
if (_tle)
|
|
||||||
*_tle = max_i + 1;
|
|
||||||
if (_gtle)
|
|
||||||
*_gtle = max_ie + 1;
|
|
||||||
if (_gscore)
|
|
||||||
*_gscore = gscore;
|
|
||||||
if (_max_off)
|
|
||||||
*_max_off = max_off;
|
|
||||||
return max;
|
|
||||||
}
|
|
||||||
|
|
||||||
typedef struct
|
|
||||||
{
|
|
||||||
int32_t h, e;
|
|
||||||
} eh_t;
|
|
||||||
|
|
||||||
int ksw_extend2_origin(int qlen, // query length 待匹配段碱基的query长度
|
|
||||||
const uint8_t *query, // read碱基序列
|
|
||||||
int tlen, // target length reference的长度
|
|
||||||
const uint8_t *target, // reference序列
|
|
||||||
int is_left, // 是不是向左扩展
|
|
||||||
int m, // 碱基种类 (5)
|
|
||||||
const int8_t *mat, // 每个位置的query和target的匹配得分 m*m
|
|
||||||
int o_del, // deletion 错配开始的惩罚系数
|
|
||||||
int e_del, // deletion extension的惩罚系数
|
|
||||||
int o_ins, // insertion 错配开始的惩罚系数
|
|
||||||
int e_ins, // insertion extension的惩罚系数
|
|
||||||
int w, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
|
||||||
int end_bonus,
|
|
||||||
int zdrop,
|
|
||||||
int h0, // 该seed的初始得分（完全匹配query的碱基数）
|
|
||||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
|
||||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
|
||||||
int *_gtle, // query全部匹配上的target的长度
|
|
||||||
int *_gscore, // query的端到端匹配得分
|
|
||||||
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
|
||||||
{
|
|
||||||
eh_t *eh; // score array
|
|
||||||
int8_t *qp; // query profile
|
|
||||||
int i, j, k, oe_del = o_del + e_del, oe_ins = o_ins + e_ins, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
|
||||||
uint8_t *qmem, *ref, *seq;
|
|
||||||
assert(h0 > 0);
|
|
||||||
// allocate memory
|
|
||||||
qp = malloc(qlen * m);
|
|
||||||
eh = calloc(qlen + 1, 8);
|
|
||||||
qmem = malloc(qlen + tlen);
|
|
||||||
seq = (uint8_t *)&qmem[0];
|
|
||||||
ref = (uint8_t *)&qmem[qlen];
|
|
||||||
if (is_left)
|
|
||||||
{
|
|
||||||
for (i = 0; i < qlen; ++i)
|
|
||||||
seq[i] = query[qlen - 1 - i];
|
|
||||||
for (i = 0; i < tlen; ++i)
|
|
||||||
ref[i] = target[tlen - 1 - i];
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
for (i = 0; i < qlen; ++i)
|
|
||||||
seq[i] = query[i];
|
|
||||||
for (i = 0; i < tlen; ++i)
|
|
||||||
ref[i] = target[i];
|
|
||||||
}
|
|
||||||
// generate the query profile
|
|
||||||
for (k = i = 0; k < m; ++k)
|
|
||||||
{
|
|
||||||
const int8_t *p = &mat[k * m];
|
|
||||||
for (j = 0; j < qlen; ++j)
|
|
||||||
qp[i++] = p[seq[j]];
|
|
||||||
}
|
|
||||||
// fill the first row
|
|
||||||
eh[0].h = h0;
|
|
||||||
eh[1].h = h0 > oe_ins ? h0 - oe_ins : 0;
|
|
||||||
for (j = 2; j <= qlen && eh[j - 1].h > e_ins; ++j)
|
|
||||||
eh[j].h = eh[j - 1].h - e_ins;
|
|
||||||
// adjust $w if it is too large
|
|
||||||
k = m * m;
|
|
||||||
for (i = 0, max = 0; i < k; ++i) // get the max score
|
|
||||||
max = max > mat[i] ? max : mat[i];
|
|
||||||
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.);
|
|
||||||
max_ins = max_ins > 1 ? max_ins : 1;
|
|
||||||
w = w < max_ins ? w : max_ins;
|
|
||||||
max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.);
|
|
||||||
max_del = max_del > 1 ? max_del : 1;
|
|
||||||
w = w < max_del ? w : max_del; // TODO: is this necessary?
|
|
||||||
// printf("%d\n", w);
|
|
||||||
// DP loop
|
|
||||||
max = h0, max_i = max_j = -1;
|
|
||||||
max_ie = -1, gscore = -1;
|
|
||||||
max_off = 0;
|
|
||||||
beg = 0, end = qlen;
|
|
||||||
|
|
||||||
for (i = 0; LIKELY(i < tlen); ++i)
|
|
||||||
{
|
|
||||||
int t, f = 0, h1, m = 0, mj = -1;
|
|
||||||
int8_t *q = &qp[ref[i] * qlen];
|
|
||||||
// apply the band and the constraint (if provided)
|
|
||||||
if (beg < i - w)
|
|
||||||
beg = i - w;
|
|
||||||
if (end > i + w + 1)
|
|
||||||
end = i + w + 1;
|
|
||||||
// if (end > qlen) end = qlen; 没用
|
|
||||||
// compute the first column
|
|
||||||
if (beg == 0)
|
|
||||||
{
|
|
||||||
h1 = h0 - (o_del + e_del * (i + 1));
|
|
||||||
if (h1 < 0)
|
|
||||||
h1 = 0;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
h1 = 0;
|
|
||||||
for (j = beg; LIKELY(j < end); ++j)
|
|
||||||
{
|
|
||||||
// At the beginning of the loop: eh[j] = { H(i-1,j-1), E(i,j) }, f = F(i,j) and h1 = H(i,j-1)
|
|
||||||
// Similar to SSE2-SW, cells are computed in the following order:
|
|
||||||
// H(i,j) = max{H(i-1,j-1)+S(i,j), E(i,j), F(i,j)}
|
|
||||||
// E(i+1,j) = max{H(i,j)-gapo, E(i,j)} - gape
|
|
||||||
// F(i,j+1) = max{H(i,j)-gapo, F(i,j)} - gape
|
|
||||||
eh_t *p = &eh[j];
|
|
||||||
int h, M = p->h, e = p->e; // get H(i-1,j-1) and E(i-1,j)
|
|
||||||
p->h = h1; // set H(i,j-1) for the next row
|
|
||||||
M = M ? M + q[j] : 0; // separating H and M to disallow a cigar like "100M3I3D20M"
|
|
||||||
h = M > e ? M : e; // e and f are guaranteed to be non-negative, so h>=0 even if M<0
|
|
||||||
h = h > f ? h : f;
|
|
||||||
h1 = h; // save H(i,j) to h1 for the next column
|
|
||||||
mj = m > h ? mj : j; // record the position where max score is achieved
|
|
||||||
m = m > h ? m : h; // m is stored at eh[mj+1]
|
|
||||||
t = M - oe_del;
|
|
||||||
t = t > 0 ? t : 0;
|
|
||||||
e -= e_del;
|
|
||||||
e = e > t ? e : t; // computed E(i+1,j)
|
|
||||||
p->e = e; // save E(i+1,j) for the next row
|
|
||||||
t = M - oe_ins;
|
|
||||||
t = t > 0 ? t : 0;
|
|
||||||
f -= e_ins;
|
|
||||||
f = f > t ? f : t; // computed F(i,j+1)
|
|
||||||
}
|
|
||||||
eh[end].h = h1;
|
|
||||||
eh[end].e = 0;
|
|
||||||
if (j == qlen)
|
|
||||||
{
|
|
||||||
max_ie = gscore > h1 ? max_ie : i;
|
|
||||||
gscore = gscore > h1 ? gscore : h1;
|
|
||||||
}
|
|
||||||
if (m == 0)
|
|
||||||
break;
|
|
||||||
if (m > max)
|
|
||||||
{
|
|
||||||
max = m, max_i = i, max_j = mj;
|
|
||||||
max_off = max_off > abs(mj - i) ? max_off : abs(mj - i);
|
|
||||||
}
|
|
||||||
else if (zdrop > 0)
|
|
||||||
{
|
|
||||||
if (i - max_i > mj - max_j)
|
|
||||||
{
|
|
||||||
if (max - m - ((i - max_i) - (mj - max_j)) * e_del > zdrop)
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
if (max - m - ((mj - max_j) - (i - max_i)) * e_ins > zdrop)
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// update beg and end for the next round
|
|
||||||
for (j = beg; LIKELY(j < end) && eh[j].h == 0 && eh[j].e == 0; ++j)
|
|
||||||
;
|
|
||||||
beg = j;
|
|
||||||
for (j = end; LIKELY(j >= beg) && eh[j].h == 0 && eh[j].e == 0; --j)
|
|
||||||
;
|
|
||||||
end = j + 2 < qlen ? j + 2 : qlen;
|
|
||||||
// beg = 0; end = qlen; // uncomment this line for debugging
|
|
||||||
}
|
|
||||||
|
|
||||||
free(eh);
|
|
||||||
free(qp);
|
|
||||||
free(qmem);
|
|
||||||
if (_qle)
|
if (_qle)
|
||||||
*_qle = max_j + 1;
|
*_qle = max_j + 1;
|
||||||
if (_tle)
|
if (_tle)
|
||||||
|
|
@ -198,28 +198,26 @@ static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = {
|
||||||
cur_match_arr = next_match_arr; \
|
cur_match_arr = next_match_arr; \
|
||||||
next_match_arr = tmp;
|
next_match_arr = tmp;
|
||||||
|
|
||||||
// uint8_t mem_addr[102400];
|
int ksw_extend_avx2_aligned(thread_mem_t *tmem,
|
||||||
|
int qlen, // query length 待匹配段碱基的query长度
|
||||||
int ksw_avx2_aligned(thread_mem_t *tmem,
|
const uint8_t *query, // read碱基序列
|
||||||
int qlen, // query length 待匹配段碱基的query长度
|
int tlen, // target length reference的长度
|
||||||
const uint8_t *query, // read碱基序列
|
const uint8_t *target, // reference序列
|
||||||
int tlen, // target length reference的长度
|
int extend_left, // 是不是向左扩展
|
||||||
const uint8_t *target, // reference序列
|
int o_del, // deletion 错配开始的惩罚系数
|
||||||
int extend_left, // 是不是向左扩展
|
int e_del, // deletion extension的惩罚系数
|
||||||
int o_del, // deletion 错配开始的惩罚系数
|
int o_ins, // insertion 错配开始的惩罚系数
|
||||||
int e_del, // deletion extension的惩罚系数
|
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||||
int o_ins, // insertion 错配开始的惩罚系数
|
int base_match_score, // 碱基match时的分数
|
||||||
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||||
int base_match_score, // 碱基match时的分数
|
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||||
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||||
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||||
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||||
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
int *_gtle, // query全部匹配上的target的长度
|
||||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
int *_gscore, // query的端到端匹配得分
|
||||||
int *_gtle, // query全部匹配上的target的长度
|
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
||||||
int *_gscore, // query的端到端匹配得分
|
|
||||||
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
|
||||||
{
|
{
|
||||||
int16_t *cur_match_arr, *next_match_arr,
|
int16_t *cur_match_arr, *next_match_arr,
|
||||||
*last_max_arr, *cur_max_arr, *next_max_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)
|
if (read_start_pos > read_end_pos)
|
||||||
break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug
|
break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug
|
||||||
|
|
||||||
// read_start_pos = 1;
|
read_start_pos = 1;
|
||||||
// read_end_pos = qlen;
|
read_end_pos = qlen;
|
||||||
ref_end_pos = di - (read_start_pos - 1); // ref开始计算的位置，倒序
|
ref_end_pos = di - (read_start_pos - 1); // ref开始计算的位置，倒序
|
||||||
span = read_end_pos - read_start_pos;
|
span = read_end_pos - read_start_pos;
|
||||||
ref_start_pos = ref_end_pos - span - 1; // 0开始的ref索引位置
|
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)
|
// for (j = read_start_pos; LIKELY(j <= read_end_pos); ++j)
|
||||||
{
|
//{
|
||||||
int has_val = cur_max_arr[j - 1] | next_max_arr[j];
|
// int has_val = cur_max_arr[j - 1] | next_max_arr[j];
|
||||||
if (has_val)
|
// if (has_val)
|
||||||
{
|
// {
|
||||||
break;
|
// break;
|
||||||
}
|
// }
|
||||||
}
|
//}
|
||||||
read_start_pos = j;
|
// read_start_pos = j;
|
||||||
|
//
|
||||||
next_max_arr[read_end_pos + 1] = 0;
|
// next_max_arr[read_end_pos + 1] = 0;
|
||||||
for (j = read_end_pos + 1; LIKELY(j >= read_start_pos); --j)
|
// for (j = read_end_pos + 1; LIKELY(j >= read_start_pos); --j)
|
||||||
{
|
//{
|
||||||
int has_val = cur_max_arr[j - 1] | next_max_arr[j];
|
// int has_val = cur_max_arr[j - 1] | next_max_arr[j];
|
||||||
if (has_val)
|
// if (has_val)
|
||||||
{
|
// {
|
||||||
break;
|
// break;
|
||||||
}
|
// }
|
||||||
}
|
//}
|
||||||
read_end_pos = j + 1 <= qlen ? j + 1 : qlen;
|
// read_end_pos = j + 1 <= qlen ? j + 1 : qlen;
|
||||||
// swap m, h, e, f
|
// swap m, h, e, f
|
||||||
SWAP_DATA_POINTER;
|
SWAP_DATA_POINTER;
|
||||||
}
|
}
|
||||||
|
|
@ -5,6 +5,7 @@
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include <immintrin.h>
|
#include <immintrin.h>
|
||||||
#include <emmintrin.h>
|
#include <emmintrin.h>
|
||||||
|
#include "thread_mem.h"
|
||||||
|
|
||||||
#ifdef __GNUC__
|
#ifdef __GNUC__
|
||||||
#define LIKELY(x) __builtin_expect((x), 1)
|
#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);
|
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 \
|
#define SIMD_STORE \
|
||||||
max_vec = _mm256_max_epu8(max_vec, hn_vec); \
|
max_vec = _mm256_max_epu8(max_vec, hn_vec); \
|
||||||
_mm256_storeu_si256((__m256i *)&next_del_arr[j], en_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; \
|
cur_match_arr = next_match_arr; \
|
||||||
next_match_arr = tmp;
|
next_match_arr = tmp;
|
||||||
|
|
||||||
uint8_t mem[102400];
|
// uint8_t mem[102400];
|
||||||
|
|
||||||
int bsw_avx2(int qlen, // query length 待匹配段碱基的query长度
|
int ksw_extend_avx2_heuristics(thread_mem_t *tmem,
|
||||||
const uint8_t *query, // read碱基序列
|
int qlen, // query length 待匹配段碱基的query长度
|
||||||
int tlen, // target length reference的长度
|
const uint8_t *query, // read碱基序列
|
||||||
const uint8_t *target, // reference序列
|
int tlen, // target length reference的长度
|
||||||
int extend_left, // 是不是向左扩展
|
const uint8_t *target, // reference序列
|
||||||
int o_del, // deletion 错配开始的惩罚系数
|
int extend_left, // 是不是向左扩展
|
||||||
int e_del, // deletion extension的惩罚系数
|
int o_del, // deletion 错配开始的惩罚系数
|
||||||
int o_ins, // insertion 错配开始的惩罚系数
|
int e_del, // deletion extension的惩罚系数
|
||||||
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
int o_ins, // insertion 错配开始的惩罚系数
|
||||||
int base_match_score, // 碱基match时的分数
|
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||||
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
int base_match_score, // 碱基match时的分数
|
||||||
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||||
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||||
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||||
int *_gtle, // query全部匹配上的target的长度
|
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||||
int *_gscore, // query的端到端匹配得分
|
int *_gtle, // query全部匹配上的target的长度
|
||||||
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
int *_gscore, // query的端到端匹配得分
|
||||||
|
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
||||||
{
|
{
|
||||||
return 0;
|
// return 0;
|
||||||
int16_t *mA, *eA, *hA, *fA,
|
int16_t *mA, *eA, *hA, *fA,
|
||||||
*cur_match_arr, *next_match_arr,
|
*cur_match_arr, *next_match_arr,
|
||||||
*last_max_arr, *cur_max_arr, *next_max_arr,
|
*last_max_arr, *cur_max_arr, *next_max_arr,
|
||||||
*cur_del_arr, *next_del_arr,
|
*cur_del_arr, *next_del_arr,
|
||||||
*cur_ins_arr, *next_ins_arr; // hA0保存上上个col的H，其他的保存上个H E F M
|
*cur_ins_arr, *next_ins_arr; // hA0保存上上个col的H，其他的保存上个H E F M
|
||||||
int16_t *seq, *ref;
|
int16_t *seq, *ref;
|
||||||
// uint8_t *mem;
|
uint8_t *mem;
|
||||||
int16_t *qtmem, *vmem;
|
int16_t *qtmem, *vmem;
|
||||||
int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH;
|
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 Dloop = tlen + qlen; // 循环跳出条件 D从1开始遍历
|
||||||
int span, beg1, end1; // 边界条件计算
|
int span, beg1, end1; // 边界条件计算
|
||||||
int col_size = qlen + 2 + SIMD_WIDTH;
|
int col_size = qlen + 2 + SIMD_WIDTH;
|
||||||
int val_mem_size = (col_size * 9 * 2 + 31) >> 5 << 5; // 32字节的整数倍
|
int val_mem_size = (col_size * 9 * 2 + 31) >> 5 << 5; // 32字节的整数倍
|
||||||
int mem_size = (seq_size + ref_size) * 2 + val_mem_size;
|
int mem_size = (seq_size + ref_size) * 2 + val_mem_size;
|
||||||
int simd_i = -1, simd_j = -1;
|
|
||||||
|
|
||||||
SIMD_INIT; // 初始化simd用的数据
|
SIMD_INIT; // 初始化simd用的数据
|
||||||
|
|
||||||
|
|
@ -245,6 +236,7 @@ int bsw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
|
|
||||||
// allocate memory
|
// allocate memory
|
||||||
// mem = malloc(mem_size);
|
// mem = malloc(mem_size);
|
||||||
|
mem = thread_mem_request(tmem, mem_size);
|
||||||
qtmem = (int16_t *)&mem[0];
|
qtmem = (int16_t *)&mem[0];
|
||||||
seq = &qtmem[0];
|
seq = &qtmem[0];
|
||||||
ref = &qtmem[seq_size];
|
ref = &qtmem[seq_size];
|
||||||
|
|
@ -313,7 +305,6 @@ int bsw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
gscore = 0;
|
gscore = 0;
|
||||||
}
|
}
|
||||||
|
|
||||||
int m_last = 0;
|
|
||||||
int iend;
|
int iend;
|
||||||
|
|
||||||
for (D = 1; LIKELY(D < Dloop); ++D)
|
for (D = 1; LIKELY(D < Dloop); ++D)
|
||||||
|
|
@ -331,8 +322,8 @@ int bsw_avx2(int qlen, // query length 待匹配段碱基的query
|
||||||
}
|
}
|
||||||
end1 = MIN(qlen, beg1 + span);
|
end1 = MIN(qlen, beg1 + span);
|
||||||
|
|
||||||
// beg = 1;
|
beg = 1;
|
||||||
// end = qlen;
|
end = qlen;
|
||||||
|
|
||||||
if (beg < beg1)
|
if (beg < beg1)
|
||||||
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;
|
max_ie = gscore > next_max_arr[qlen] ? max_ie : iStart;
|
||||||
gscore = gscore > next_max_arr[qlen] ? gscore : next_max_arr[qlen];
|
gscore = gscore > next_max_arr[qlen] ? gscore : next_max_arr[qlen];
|
||||||
}
|
}
|
||||||
// if (m == 0 && m_last == 0)
|
|
||||||
// break; // 一定要注意，斜对角遍历和按列遍历的不同点
|
|
||||||
if (m > max)
|
if (m > max)
|
||||||
{
|
{
|
||||||
max = m, max_i = mi, max_j = mj;
|
max = m, max_i = mi, max_j = mj;
|
||||||
max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi);
|
max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi);
|
||||||
}
|
}
|
||||||
|
|
||||||
// 调整计算的边界
|
// for (j = beg; LIKELY(j <= end); ++j)
|
||||||
|
|
||||||
// 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)
|
|
||||||
//{
|
//{
|
||||||
// __m256i h1 = _mm256_loadu_si256((__m256i *)(&cur_max_arr[j - 1]));
|
// int has_val = cur_max_arr[j - 1] | next_max_arr[j];
|
||||||
// __m256i h2 = _mm256_loadu_si256((__m256i *)(&next_max_arr[j]));
|
// if (has_val)
|
||||||
// __m256i orvec = _mm256_or_si256(h1, h2);
|
// {
|
||||||
// int *val = (int *)&orvec;
|
// break;
|
||||||
// for (i = 0; i < SIMD_WIDTH; ++i)
|
// }
|
||||||
// if (val[i])
|
// }
|
||||||
// {
|
|
||||||
// pos = SIMD_WIDTH - 1 - i;
|
|
||||||
// break;
|
|
||||||
// }
|
|
||||||
//}
|
|
||||||
// beg = j;
|
// beg = j;
|
||||||
for (j = beg; LIKELY(j <= end); ++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)
|
// int has_val = cur_max_arr[j - 1] | next_max_arr[j];
|
||||||
{
|
// if (has_val)
|
||||||
break;
|
// {
|
||||||
}
|
// break;
|
||||||
}
|
// }
|
||||||
beg = j;
|
// }
|
||||||
|
// 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 m, h, e, f
|
||||||
SWAP_DATA_POINTER;
|
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);
|
// free(mem);
|
||||||
|
thread_mem_release(tmem, mem_size);
|
||||||
if (_qle)
|
if (_qle)
|
||||||
*_qle = max_j + 1;
|
*_qle = max_j + 1;
|
||||||
if (_tle)
|
if (_tle)
|
||||||
|
|
@ -0,0 +1,434 @@
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <assert.h>
|
||||||
|
#include <emmintrin.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <immintrin.h>
|
||||||
|
#include <emmintrin.h>
|
||||||
|
#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;
|
||||||
|
}
|
||||||
|
|
@ -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, 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}};
|
{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};
|
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 _MM_SHUFFLE(0, 1, 2, 3)
|
||||||
#define permute_mask 27
|
#define permute_mask 27
|
||||||
// 初始化变量
|
// 初始化变量
|
||||||
#define SIMD_INIT \
|
#define SIMD_INIT \
|
||||||
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
||||||
__m256i zero_vec; \
|
__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_del_vec; \
|
||||||
__m256i oe_ins_vec; \
|
__m256i oe_ins_vec; \
|
||||||
__m256i e_del_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); \
|
oe_ins_vec = _mm256_set1_epi8(oe_ins); \
|
||||||
e_del_vec = _mm256_set1_epi8(e_del); \
|
e_del_vec = _mm256_set1_epi8(e_del); \
|
||||||
e_ins_vec = _mm256_set1_epi8(e_ins); \
|
e_ins_vec = _mm256_set1_epi8(e_ins); \
|
||||||
__m256i match_sc_vec = _mm256_set1_epi8(a); \
|
__m256i match_sc_vec = _mm256_set1_epi8(base_match_score); \
|
||||||
__m256i mis_sc_vec = _mm256_set1_epi8(b); \
|
__m256i mis_sc_vec = _mm256_set1_epi8(base_mis_score); \
|
||||||
__m256i amb_sc_vec = _mm256_set1_epi8(1); \
|
__m256i amb_sc_vec = _mm256_set1_epi8(1); \
|
||||||
__m256i amb_vec = _mm256_set1_epi8(4); \
|
__m256i amb_vec = _mm256_set1_epi8(4); \
|
||||||
reverse_mask_vec = _mm256_loadu_si256((__m256i *)(reverse_mask)); \
|
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]) \
|
if (read_seq[mj + 2] == ref_seq[mi + 1 + SIMD_WIDTH]) \
|
||||||
{ \
|
{ \
|
||||||
m += a; \
|
m += base_match_score; \
|
||||||
} \
|
} \
|
||||||
else \
|
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; \
|
mA1 = mA2; \
|
||||||
mA2 = tmp;
|
mA2 = tmp;
|
||||||
|
|
||||||
int ksw_avx2_u8_aligned(thread_mem_t *tmem,
|
int ksw_extend_avx2_u8_aligned(thread_mem_t *tmem,
|
||||||
int qlen, // query length 待匹配段碱基的query长度
|
int qlen, // query length 待匹配段碱基的query长度
|
||||||
const uint8_t *query, // read碱基序列
|
const uint8_t *query, // read碱基序列
|
||||||
int tlen, // target length reference的长度
|
int tlen, // target length reference的长度
|
||||||
const uint8_t *target, // reference序列
|
const uint8_t *target, // reference序列
|
||||||
int is_left, // 是不是向左扩展
|
int extend_left, // 是不是向左扩展
|
||||||
int m, // 碱基种类 (5)
|
int o_del, // deletion 错配开始的惩罚系数
|
||||||
const int8_t *mat, // 每个位置的query和target的匹配得分 m*m
|
int e_del, // deletion extension的惩罚系数
|
||||||
int o_del, // deletion 错配开始的惩罚系数
|
int o_ins, // insertion 错配开始的惩罚系数
|
||||||
int e_del, // deletion extension的惩罚系数
|
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||||
int o_ins, // insertion 错配开始的惩罚系数
|
int base_match_score, // 碱基match时的分数
|
||||||
int e_ins, // insertion extension的惩罚系数
|
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||||
int a, // 碱基match时的分数
|
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||||
int b, // 碱基mismatch时的惩罚分数（正数）
|
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||||
int w, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||||
int end_bonus,
|
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||||
int zdrop,
|
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||||
int h0, // 该seed的初始得分（完全匹配query的碱基数）
|
int *_gtle, // query全部匹配上的target的长度
|
||||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
int *_gscore, // query的端到端匹配得分
|
||||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
int *_max_off) // 取得最大得分时在query和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;
|
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 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;
|
uint8_t *mem_addr;
|
||||||
int read_size = align_number(qlen * BASE_BYTES + MEM_ALIGN_BYTES);
|
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用的数据
|
SIMD_INIT; // 初始化simd用的数据
|
||||||
|
|
||||||
assert(h0 > 0);
|
assert(init_score > 0);
|
||||||
|
|
||||||
mem_addr = thread_mem_request(tmem, request_mem_size);
|
mem_addr = thread_mem_request(tmem, request_mem_size);
|
||||||
mem_addr = (void *)align_mem((uint64_t)mem_addr);
|
mem_addr = (void *)align_mem((uint64_t)mem_addr);
|
||||||
ref_seq = (uint8_t *)&mem_addr[0];
|
ref_seq = (uint8_t *)&mem_addr[0];
|
||||||
read_seq = (uint8_t *)(mem_addr + ref_size + SIMD_BYTES - BASE_BYTES);
|
read_seq = (uint8_t *)(mem_addr + ref_size + SIMD_BYTES - BASE_BYTES);
|
||||||
|
|
||||||
if (is_left)
|
if (extend_left)
|
||||||
{
|
{
|
||||||
for (i = 0; i < qlen; ++i)
|
for (i = 0; i < qlen; ++i)
|
||||||
read_seq[i + 1] = query[qlen - 1 - 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];
|
fA1 = &eA2[score_element_num];
|
||||||
fA2 = &fA1[score_element_num];
|
fA2 = &fA1[score_element_num];
|
||||||
|
|
||||||
// adjust $w if it is too large
|
// adjust $window_size if it is too large
|
||||||
k = m * m;
|
|
||||||
// get the max score
|
max = base_match_score;
|
||||||
for (i = 0, max = 0; i < k; ++i)
|
|
||||||
max = max > mat[i] ? max : mat[i];
|
|
||||||
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.);
|
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.);
|
||||||
max_ins = max_ins > 1 ? max_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 = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.);
|
||||||
max_del = max_del > 1 ? max_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)
|
if (tlen < qlen)
|
||||||
w = MIN(tlen - 1, w);
|
window_size = MIN(tlen - 1, window_size);
|
||||||
|
|
||||||
// DP loop
|
// DP loop
|
||||||
max = h0, max_i = max_j = -1;
|
max = init_score, max_i = max_j = -1;
|
||||||
max_ie = -1, gscore = -1;
|
max_ie = -1, gscore = -1;
|
||||||
;
|
;
|
||||||
max_off = 0;
|
max_off = 0;
|
||||||
beg = 1;
|
beg = 1;
|
||||||
end = qlen;
|
end = qlen;
|
||||||
// init h0
|
// init init_score
|
||||||
hA0[0] = h0; // 左上角
|
hA0[0] = init_score; // 左上角
|
||||||
|
|
||||||
if (qlen == 0 || tlen == 0)
|
if (qlen == 0 || tlen == 0)
|
||||||
back_diagnal_num = 0; // 防止意外情况
|
back_diagnal_num = 0; // 防止意外情况
|
||||||
if (w >= qlen)
|
if (window_size >= qlen)
|
||||||
{
|
{
|
||||||
max_ie = 0;
|
max_ie = 0;
|
||||||
gscore = 0;
|
gscore = 0;
|
||||||
}
|
}
|
||||||
|
|
||||||
int m_last = 0;
|
|
||||||
int iend;
|
int iend;
|
||||||
|
|
||||||
for (D = 1; LIKELY(D < back_diagnal_num); ++D)
|
for (D = 1; LIKELY(D < back_diagnal_num); ++D)
|
||||||
|
|
@ -350,13 +339,13 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem,
|
||||||
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
||||||
if (D > tlen)
|
if (D > tlen)
|
||||||
{
|
{
|
||||||
span = MIN(back_diagnal_num - D, w);
|
span = MIN(back_diagnal_num - D, window_size);
|
||||||
beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1);
|
beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1);
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
span = MIN(D - 1, w);
|
span = MIN(D - 1, window_size);
|
||||||
beg1 = MAX(1, ((D - w) / 2) + 1);
|
beg1 = MAX(1, ((D - window_size) / 2) + 1);
|
||||||
}
|
}
|
||||||
end1 = MIN(qlen, beg1 + span);
|
end1 = MIN(qlen, beg1 + span);
|
||||||
|
|
||||||
|
|
@ -367,8 +356,8 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem,
|
||||||
if (beg > end)
|
if (beg > end)
|
||||||
break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug
|
break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug
|
||||||
|
|
||||||
// beg = 1;
|
beg = 1;
|
||||||
// end = qlen;
|
end = qlen;
|
||||||
iend = D - (beg - 1); // ref开始计算的位置，倒序
|
iend = D - (beg - 1); // ref开始计算的位置，倒序
|
||||||
span = end - beg;
|
span = end - beg;
|
||||||
iStart = iend - span - 1; // 0开始的ref索引位置
|
iStart = iend - span - 1; // 0开始的ref索引位置
|
||||||
|
|
@ -381,12 +370,12 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem,
|
||||||
// 左边界 处理f (insert)
|
// 左边界 处理f (insert)
|
||||||
if (iStart == 0)
|
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)
|
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
|
else
|
||||||
{
|
{
|
||||||
|
|
@ -438,26 +427,22 @@ int ksw_avx2_u8_aligned(thread_mem_t *tmem,
|
||||||
}
|
}
|
||||||
|
|
||||||
// 调整计算的边界
|
// 调整计算的边界
|
||||||
for (j = beg; LIKELY(j <= end); ++j)
|
// for (j = beg; LIKELY(j <= end); ++j)
|
||||||
{
|
//{
|
||||||
int has_val = hA1[j - 1] | hA2[j];
|
// int has_val = hA1[j - 1] | hA2[j];
|
||||||
if (has_val)
|
// if (has_val)
|
||||||
break;
|
// break;
|
||||||
}
|
//}
|
||||||
beg = j;
|
// beg = j;
|
||||||
hA2[end + 1] = 0;
|
// hA2[end + 1] = 0;
|
||||||
for (j = end + 1; LIKELY(j >= beg); --j)
|
// for (j = end + 1; LIKELY(j >= beg); --j)
|
||||||
{
|
//{
|
||||||
int has_val = hA1[j - 1] | hA2[j];
|
// int has_val = hA1[j - 1] | hA2[j];
|
||||||
if (has_val)
|
// if (has_val)
|
||||||
break;
|
// break;
|
||||||
}
|
//}
|
||||||
end = j + 1 <= qlen ? j + 1 : qlen;
|
// end = j + 1 <= qlen ? j + 1 : qlen;
|
||||||
|
|
||||||
// beg = 0;
|
|
||||||
// end = qlen;
|
|
||||||
|
|
||||||
m_last = m;
|
|
||||||
// swap m, h, e, f
|
// swap m, h, e, f
|
||||||
SWAP_DATA_POINTER;
|
SWAP_DATA_POINTER;
|
||||||
}
|
}
|
||||||
|
|
@ -5,6 +5,7 @@
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include <immintrin.h>
|
#include <immintrin.h>
|
||||||
#include <emmintrin.h>
|
#include <emmintrin.h>
|
||||||
|
#include "thread_mem.h"
|
||||||
|
|
||||||
#ifdef __GNUC__
|
#ifdef __GNUC__
|
||||||
#define LIKELY(x) __builtin_expect((x), 1)
|
#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 \
|
#define SIMD_INIT \
|
||||||
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
||||||
__m256i zero_vec; \
|
__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_del_vec; \
|
||||||
__m256i oe_ins_vec; \
|
__m256i oe_ins_vec; \
|
||||||
__m256i e_del_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); \
|
oe_ins_vec = _mm256_set1_epi8(oe_ins); \
|
||||||
e_del_vec = _mm256_set1_epi8(e_del); \
|
e_del_vec = _mm256_set1_epi8(e_del); \
|
||||||
e_ins_vec = _mm256_set1_epi8(e_ins); \
|
e_ins_vec = _mm256_set1_epi8(e_ins); \
|
||||||
__m256i match_sc_vec = _mm256_set1_epi8(a); \
|
__m256i match_sc_vec = _mm256_set1_epi8(base_match_score); \
|
||||||
__m256i mis_sc_vec = _mm256_set1_epi8(b); \
|
__m256i mis_sc_vec = _mm256_set1_epi8(base_mis_score); \
|
||||||
__m256i amb_sc_vec = _mm256_set1_epi8(1); \
|
__m256i amb_sc_vec = _mm256_set1_epi8(1); \
|
||||||
__m256i amb_vec = _mm256_set1_epi8(4); \
|
__m256i amb_vec = _mm256_set1_epi8(4); \
|
||||||
reverse_mask_vec = _mm256_loadu_si256((__m256i *)(reverse_mask)); \
|
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]) \
|
if (seq[mj + 1] == ref[mi + 1 + SIMD_WIDTH]) \
|
||||||
{ \
|
{ \
|
||||||
m += a; \
|
m += base_match_score; \
|
||||||
} \
|
} \
|
||||||
else \
|
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; \
|
mA1 = mA2; \
|
||||||
mA2 = tmp;
|
mA2 = tmp;
|
||||||
|
|
||||||
uint8_t mem1[102400];
|
// uint8_t mem1[102400];
|
||||||
|
|
||||||
int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的query长度
|
int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem,
|
||||||
const uint8_t *query, // read碱基序列
|
int qlen, // query length 待匹配段碱基的query长度
|
||||||
int tlen, // target length reference的长度
|
const uint8_t *query, // read碱基序列
|
||||||
const uint8_t *target, // reference序列
|
int tlen, // target length reference的长度
|
||||||
int is_left, // 是不是向左扩展
|
const uint8_t *target, // reference序列
|
||||||
int m, // 碱基种类 (5)
|
int extend_left, // 是不是向左扩展
|
||||||
const int8_t *mat, // 每个位置的query和target的匹配得分 m*m
|
int o_del, // deletion 错配开始的惩罚系数
|
||||||
int o_del, // deletion 错配开始的惩罚系数
|
int e_del, // deletion extension的惩罚系数
|
||||||
int e_del, // deletion extension的惩罚系数
|
int o_ins, // insertion 错配开始的惩罚系数
|
||||||
int o_ins, // insertion 错配开始的惩罚系数
|
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||||
int e_ins, // insertion extension的惩罚系数
|
int base_match_score, // 碱基match时的分数
|
||||||
int a, // 碱基match时的分数
|
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||||
int b, // 碱基mismatch时的惩罚分数（正数）
|
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||||
int w, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||||
int end_bonus,
|
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||||
int zdrop,
|
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||||
int h0, // 该seed的初始得分（完全匹配query的碱基数）
|
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
int *_gtle, // query全部匹配上的target的长度
|
||||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
int *_gscore, // query的端到端匹配得分
|
||||||
int *_gtle, // query全部匹配上的target的长度
|
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
||||||
int *_gscore, // query的端到端匹配得分
|
|
||||||
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
|
||||||
{
|
{
|
||||||
// return h0;
|
// return init_score;
|
||||||
#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
|
|
||||||
uint8_t *mA, *hA, *eA, *fA, *mA1, *mA2, *hA0, *hA1, *eA1, *fA1, *hA2, *eA2, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M
|
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 *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 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 Dloop = tlen + qlen; // 循环跳出条件
|
||||||
int span, beg1, end1; // 边界条件计算
|
int span, beg1, end1; // 边界条件计算
|
||||||
int col_size = qlen + 2 + SIMD_WIDTH;
|
int col_size = qlen + 2 + SIMD_WIDTH;
|
||||||
|
|
@ -260,14 +251,15 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que
|
||||||
|
|
||||||
SIMD_INIT; // 初始化simd用的数据
|
SIMD_INIT; // 初始化simd用的数据
|
||||||
|
|
||||||
assert(h0 > 0);
|
assert(init_score > 0);
|
||||||
|
|
||||||
// allocate memory
|
// allocate memory
|
||||||
// mem = malloc(mem_size);
|
// mem = malloc(mem_size);
|
||||||
|
mem = thread_mem_request(tmem, mem_size);
|
||||||
qtmem = &mem[0];
|
qtmem = &mem[0];
|
||||||
seq = (uint8_t *)&qtmem[0];
|
seq = (uint8_t *)&qtmem[0];
|
||||||
ref = (uint8_t *)&qtmem[seq_size];
|
ref = (uint8_t *)&qtmem[seq_size];
|
||||||
if (is_left)
|
if (extend_left)
|
||||||
{
|
{
|
||||||
for (i = 0; i < qlen; ++i)
|
for (i = 0; i < qlen; ++i)
|
||||||
seq[i] = query[qlen - 1 - i];
|
seq[i] = query[qlen - 1 - i];
|
||||||
|
|
@ -303,39 +295,36 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que
|
||||||
fA1 = &fA[0];
|
fA1 = &fA[0];
|
||||||
fA2 = &fA[col_size];
|
fA2 = &fA[col_size];
|
||||||
|
|
||||||
// adjust $w if it is too large
|
// adjust $window_size if it is too large
|
||||||
k = m * m;
|
|
||||||
// get the max score
|
// get the max score
|
||||||
for (i = 0, max = 0; i < k; ++i)
|
max = base_match_score;
|
||||||
max = max > mat[i] ? max : mat[i];
|
|
||||||
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.);
|
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.);
|
||||||
max_ins = max_ins > 1 ? max_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 = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.);
|
||||||
max_del = max_del > 1 ? max_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)
|
if (tlen < qlen)
|
||||||
w = MIN(tlen - 1, w);
|
window_size = MIN(tlen - 1, window_size);
|
||||||
|
|
||||||
// DP loop
|
// DP loop
|
||||||
max = h0, max_i = max_j = -1;
|
max = init_score, max_i = max_j = -1;
|
||||||
max_ie = -1, gscore = -1;
|
max_ie = -1, gscore = -1;
|
||||||
;
|
;
|
||||||
max_off = 0;
|
max_off = 0;
|
||||||
beg = 1;
|
beg = 1;
|
||||||
end = qlen;
|
end = qlen;
|
||||||
// init h0
|
// init init_score
|
||||||
hA0[0] = h0; // 左上角
|
hA0[0] = init_score; // 左上角
|
||||||
|
|
||||||
if (qlen == 0 || tlen == 0)
|
if (qlen == 0 || tlen == 0)
|
||||||
Dloop = 0; // 防止意外情况
|
Dloop = 0; // 防止意外情况
|
||||||
if (w >= qlen)
|
if (window_size >= qlen)
|
||||||
{
|
{
|
||||||
max_ie = 0;
|
max_ie = 0;
|
||||||
gscore = 0;
|
gscore = 0;
|
||||||
}
|
}
|
||||||
|
|
||||||
int m_last = 0;
|
|
||||||
int iend;
|
int iend;
|
||||||
#ifdef SHOW_PERF
|
#ifdef SHOW_PERF
|
||||||
// time_bsw_init += get_mseconds() - start_time;
|
// time_bsw_init += get_mseconds() - start_time;
|
||||||
|
|
@ -349,13 +338,13 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que
|
||||||
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
||||||
if (D > tlen)
|
if (D > tlen)
|
||||||
{
|
{
|
||||||
span = MIN(Dloop - D, w);
|
span = MIN(Dloop - D, window_size);
|
||||||
beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1);
|
beg1 = MAX(D - tlen + 1, ((D - window_size) / 2) + 1);
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
span = MIN(D - 1, w);
|
span = MIN(D - 1, window_size);
|
||||||
beg1 = MAX(1, ((D - w) / 2) + 1);
|
beg1 = MAX(1, ((D - window_size) / 2) + 1);
|
||||||
}
|
}
|
||||||
end1 = MIN(qlen, beg1 + span);
|
end1 = MIN(qlen, beg1 + span);
|
||||||
|
|
||||||
|
|
@ -366,8 +355,8 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que
|
||||||
if (beg > end)
|
if (beg > end)
|
||||||
break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug
|
break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug
|
||||||
|
|
||||||
// beg = 1;
|
beg = 1;
|
||||||
// end = qlen;
|
end = qlen;
|
||||||
iend = D - (beg - 1); // ref开始计算的位置，倒序
|
iend = D - (beg - 1); // ref开始计算的位置，倒序
|
||||||
span = end - beg;
|
span = end - beg;
|
||||||
iStart = iend - span - 1; // 0开始的ref索引位置
|
iStart = iend - span - 1; // 0开始的ref索引位置
|
||||||
|
|
@ -380,12 +369,12 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que
|
||||||
// 左边界 处理f (insert)
|
// 左边界 处理f (insert)
|
||||||
if (iStart == 0)
|
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)
|
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
|
else
|
||||||
{
|
{
|
||||||
|
|
@ -418,19 +407,9 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que
|
||||||
// 存储结果
|
// 存储结果
|
||||||
SIMD_STORE;
|
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;
|
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的值
|
||||||
j = end + 1;
|
j = end + 1;
|
||||||
|
|
||||||
|
|
@ -439,50 +418,30 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que
|
||||||
max_ie = gscore > hA2[qlen] ? max_ie : iStart;
|
max_ie = gscore > hA2[qlen] ? max_ie : iStart;
|
||||||
gscore = gscore > hA2[qlen] ? gscore : hA2[qlen];
|
gscore = gscore > hA2[qlen] ? gscore : hA2[qlen];
|
||||||
}
|
}
|
||||||
// if (m == 0 && m_last == 0)
|
|
||||||
// break; // 一定要注意，斜对角遍历和按列遍历的不同点
|
|
||||||
if (m > max)
|
if (m > max)
|
||||||
{
|
{
|
||||||
max = m, max_i = mi, max_j = mj;
|
max = m, max_i = mi, max_j = mj;
|
||||||
max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi);
|
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)
|
// for (j = beg; LIKELY(j <= end); ++j)
|
||||||
{
|
//{
|
||||||
int has_val = hA1[j - 1] | hA2[j];
|
// int has_val = hA1[j - 1] | hA2[j];
|
||||||
if (has_val)
|
// if (has_val)
|
||||||
break;
|
// break;
|
||||||
}
|
//}
|
||||||
beg = j;
|
// beg = j;
|
||||||
hA2[end + 1] = 0;
|
// hA2[end + 1] = 0;
|
||||||
for (j = end + 1; LIKELY(j >= beg); --j)
|
// for (j = end + 1; LIKELY(j >= beg); --j)
|
||||||
{
|
//{
|
||||||
int has_val = hA1[j - 1] | hA2[j];
|
// int has_val = hA1[j - 1] | hA2[j];
|
||||||
if (has_val)
|
// if (has_val)
|
||||||
break;
|
// break;
|
||||||
// else
|
//}
|
||||||
// hA0[j - 1] = 0;
|
// end = j + 1 <= qlen ? j + 1 : qlen;
|
||||||
}
|
|
||||||
end = j + 1 <= qlen ? j + 1 : qlen;
|
|
||||||
|
|
||||||
// beg = 0;
|
|
||||||
// end = qlen;
|
|
||||||
|
|
||||||
m_last = m;
|
|
||||||
// swap m, h, e, f
|
// swap m, h, e, f
|
||||||
SWAP_DATA_POINTER;
|
SWAP_DATA_POINTER;
|
||||||
#ifdef SHOW_PERF
|
#ifdef SHOW_PERF
|
||||||
|
|
@ -493,6 +452,7 @@ int ksw_avx2_u8(int qlen, // query length 待匹配段碱基的que
|
||||||
}
|
}
|
||||||
|
|
||||||
// free(mem);
|
// free(mem);
|
||||||
|
thread_mem_release(tmem, mem_size);
|
||||||
if (_qle)
|
if (_qle)
|
||||||
*_qle = max_j + 1;
|
*_qle = max_j + 1;
|
||||||
if (_tle)
|
if (_tle)
|
||||||
|
|
@ -2,6 +2,7 @@
|
||||||
#include <stdlib.h>
|
#include <stdlib.h>
|
||||||
#include <assert.h>
|
#include <assert.h>
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
|
#include "thread_mem.h"
|
||||||
|
|
||||||
#ifdef __GNUC__
|
#ifdef __GNUC__
|
||||||
#define LIKELY(x) __builtin_expect((x), 1)
|
#define LIKELY(x) __builtin_expect((x), 1)
|
||||||
|
|
@ -16,15 +17,17 @@ typedef struct
|
||||||
int32_t h, e;
|
int32_t h, e;
|
||||||
} eh_t;
|
} 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
|
eh_t *eh; // score array
|
||||||
int8_t *qp; // query profile
|
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;
|
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);
|
assert(h0 > 0);
|
||||||
qp = malloc(qlen * m);
|
// qp = malloc(qlen * m);
|
||||||
eh = calloc(qlen + 1, 8);
|
// 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
|
// generate the query profile
|
||||||
for (k = i = 0; k < m; ++k)
|
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, "\n");
|
||||||
// fprintf(stderr, "%d\n", end);
|
// fprintf(stderr, "%d\n", end);
|
||||||
}
|
}
|
||||||
free(eh);
|
// free(eh);
|
||||||
free(qp);
|
// free(qp);
|
||||||
|
thread_mem_release(tmem, qlen * m + (qlen + 1) * 8);
|
||||||
if (_qle)
|
if (_qle)
|
||||||
*_qle = max_j + 1;
|
*_qle = max_j + 1;
|
||||||
if (_tle)
|
if (_tle)
|
||||||
419
main.c
419
main.c
|
|
@ -6,6 +6,8 @@
|
||||||
#include <time.h>
|
#include <time.h>
|
||||||
#include "sys/time.h"
|
#include "sys/time.h"
|
||||||
#include "thread_mem.h"
|
#include "thread_mem.h"
|
||||||
|
#include "ksw_ext.h"
|
||||||
|
#include "utils.h"
|
||||||
|
|
||||||
#define SW_NORMAL 0
|
#define SW_NORMAL 0
|
||||||
#define SW_AVX2 1
|
#define SW_AVX2 1
|
||||||
|
|
@ -15,9 +17,12 @@
|
||||||
#define BLOCK_BUF_SIZE 1048576
|
#define BLOCK_BUF_SIZE 1048576
|
||||||
#define READ_BUF_SIZE 2048
|
#define READ_BUF_SIZE 2048
|
||||||
#define SEQ_BUF_SIZE (BLOCK_BUF_SIZE + READ_BUF_SIZE)
|
#define SEQ_BUF_SIZE (BLOCK_BUF_SIZE + READ_BUF_SIZE)
|
||||||
|
#define INIT_ALLOC_SIZE 4096
|
||||||
|
|
||||||
#define DIVIDE_BY (CLOCKS_PER_SEC * 1.0)
|
#define DIVIDE_BY (CLOCKS_PER_SEC * 1.0)
|
||||||
|
|
||||||
|
#define KERNEL_NUM 7
|
||||||
|
|
||||||
#ifdef SHOW_PERF
|
#ifdef SHOW_PERF
|
||||||
// 用来调试，计算感兴趣部分的运行时间
|
// 用来调试，计算感兴趣部分的运行时间
|
||||||
// 获取当前毫秒数
|
// 获取当前毫秒数
|
||||||
|
|
@ -29,60 +34,89 @@ int64_t get_mseconds()
|
||||||
return clock();
|
return clock();
|
||||||
}
|
}
|
||||||
|
|
||||||
int64_t time_sw_normal = 0,
|
int64_t time_sw[KERNEL_NUM] = {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;
|
|
||||||
|
|
||||||
#endif
|
#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);
|
#define _PERFORMANCE_TEST_NORMAL(kernel_num, func) \
|
||||||
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,
|
cur_query_pos = 0; \
|
||||||
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);
|
cur_target_pos = 0; \
|
||||||
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,
|
for (i = 0; i < block_line_num; ++i) \
|
||||||
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);
|
{ \
|
||||||
|
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] = {
|
#define _PERFORMANCE_TEST_AVX2(kernel_num, func) \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
cur_query_pos = 0; \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
cur_target_pos = 0; \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5 /*'-'*/, 4, 4,
|
for (i = 0; i < block_line_num; ++i) \
|
||||||
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,
|
score[kernel_num] = func( \
|
||||||
4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
&tmem[kernel_num], \
|
||||||
4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
|
info_arr[i][0], \
|
||||||
4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
(uint8_t *)query_arr + cur_query_pos, \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
info_arr[i][1], \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
(uint8_t *)target_arr + cur_target_pos, \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
0, 6, 1, 6, 1, \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
1, 4, \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
100, 5, \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
info_arr[i][2], \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
|
&qle, &tle, >le, &gscore, &max_off[0]); \
|
||||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
|
score_total[kernel_num] += score[kernel_num]; \
|
||||||
char t_2bit2char[5] = {'A', 'C', 'G', 'T'};
|
cur_query_pos += info_arr[i][0]; \
|
||||||
// 将碱基字符转成2位编码
|
cur_target_pos += info_arr[i][1]; \
|
||||||
void convert_char_to_2bit(char *str)
|
}
|
||||||
{
|
|
||||||
int i;
|
|
||||||
for (i = 0; i < strlen(str); ++i)
|
|
||||||
str[i] = nst_nt4_table[str[i]];
|
|
||||||
}
|
|
||||||
|
|
||||||
// 读取测试数据
|
#ifdef SHOW_PERF
|
||||||
int read_data()
|
#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[])
|
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,
|
int8_t mat[25] = {1, -4, -4, -4, -1,
|
||||||
-4, 1, -4, -4, -1,
|
-4, 1, -4, -4, -1,
|
||||||
|
|
@ -91,52 +125,41 @@ int main(int argc, char *argv[])
|
||||||
-1, -1, -1, -1, -1};
|
-1, -1, -1, -1, -1};
|
||||||
int max_off[2];
|
int max_off[2];
|
||||||
int qle, tle, gtle, gscore;
|
int qle, tle, gtle, gscore;
|
||||||
thread_mem_t tmem, tmem_u8;
|
thread_mem_t tmem[KERNEL_NUM];
|
||||||
init_thread_mem(&tmem);
|
int i, j;
|
||||||
init_thread_mem(&tmem_u8);
|
for (i = 0; i < KERNEL_NUM; ++i)
|
||||||
// thread_mem_init_alloc(&tmem_u8, 10960);
|
{
|
||||||
|
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 *query_arr = (char *)malloc(SEQ_BUF_SIZE);
|
||||||
char *target_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_buf = (int *)malloc(SEQ_BUF_SIZE * sizeof(int));
|
||||||
int **info_arr = (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;
|
FILE *query_f = 0, *target_f = 0, *info_f = 0;
|
||||||
// const char *qf_path = "q.fa";
|
FILE *normal_out_f = 0, *avx2_out_f = 0, *avx2_u8_out_f = 0;
|
||||||
// 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";
|
|
||||||
query_f = fopen(qf_path, "r");
|
query_f = fopen(qf_path, "r");
|
||||||
target_f = fopen(tf_path, "r");
|
target_f = fopen(tf_path, "r");
|
||||||
info_f = fopen(if_path, "r");
|
info_f = fopen(if_path, "r");
|
||||||
|
|
||||||
// 将每次比对的得分等信息写入文件，进行debug
|
// 将每次比对的得分等信息写入文件，进行debug
|
||||||
FILE *normal_out_f = fopen("normal_out.txt", "w");
|
// normal_out_f = fopen("normal_out.txt", "w");
|
||||||
FILE *avx2_out_f = fopen("avx2_out.txt", "w");
|
// avx2_out_f = fopen("avx2_out.txt", "w");
|
||||||
FILE *avx2_u8_out_f = fopen("avx2_u8_out.txt", "w");
|
// avx2_u8_out_f = fopen("avx2_u8_out.txt", "w");
|
||||||
FILE *bsw_avx2_out_f = fopen("bsw_avx2_out.txt", "w");
|
|
||||||
|
|
||||||
// 每次读取一定量的数据，然后执行，直到处理完所有数据
|
// 每次读取一定量的数据，然后执行，直到处理完所有数据
|
||||||
int total_line_num = 0; // 目前处理的总的数据行数
|
int total_line_num = 0; // 目前处理的总的数据行数
|
||||||
int block_line_num = 0; // 当前循环包含的数据行数
|
int block_line_num = 0; // 当前循环包含的数据行数
|
||||||
int i, j;
|
int cur_query_pos, cur_target_pos;
|
||||||
// const int max_read = READ_BUF_SIZE; // 每次最多读取的字符
|
int64_t start_time;
|
||||||
char read_buf[READ_BUF_SIZE]; // 读文件缓存
|
char read_buf[READ_BUF_SIZE]; // 读文件缓存
|
||||||
// int ret_code = 0;
|
|
||||||
// 初始化info_arr数组
|
// 初始化info_arr数组
|
||||||
i = 0;
|
i = 0;
|
||||||
j = 0;
|
j = 0;
|
||||||
|
|
@ -147,55 +170,30 @@ int main(int argc, char *argv[])
|
||||||
info_arr[i] = &info_buf[j];
|
info_arr[i] = &info_buf[j];
|
||||||
i += 1;
|
i += 1;
|
||||||
j += 3;
|
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))
|
while (!feof(target_f))
|
||||||
{
|
{
|
||||||
// fprintf(stderr, "debug\n");
|
block_line_num = 0; // 记录每次读取的行数
|
||||||
block_line_num = 0;
|
|
||||||
// target序列一般占用存储最多，先读取target，看一个buf能读多少行，query和info就按照这个行数来读
|
// target序列一般占用存储最多，先读取target，看一个buf能读多少行，query和info就按照这个行数来读
|
||||||
int cur_read_size = 0;
|
int cur_read_size = 0;
|
||||||
while (!feof(target_f) && cur_read_size < BLOCK_BUF_SIZE)
|
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;
|
break;
|
||||||
int line_size = strlen(read_buf);
|
cur_read_size += line_size;
|
||||||
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);
|
|
||||||
++block_line_num;
|
++block_line_num;
|
||||||
++total_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
|
// 读query
|
||||||
cur_read_size = 0;
|
cur_read_size = 0;
|
||||||
for (i = 0; i < block_line_num; ++i)
|
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;
|
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;
|
cur_read_size += line_size;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
@ -215,199 +213,42 @@ int main(int argc, char *argv[])
|
||||||
|
|
||||||
// 性能测试
|
// 性能测试
|
||||||
|
|
||||||
int cur_query_pos = 0;
|
// normal sw
|
||||||
int cur_target_pos = 0;
|
PERFORMANCE_TEST_NORMAL(0, ksw_extend_normal);
|
||||||
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]);
|
|
||||||
|
|
||||||
// #ifdef SHOW_PERF
|
// avx2
|
||||||
// start_time = get_mseconds();
|
PERFORMANCE_TEST_AVX2(1, ksw_extend_avx2);
|
||||||
// #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]);
|
|
||||||
|
|
||||||
#ifdef SHOW_PERF
|
// avx2 u8
|
||||||
start_time = get_mseconds();
|
PERFORMANCE_TEST_AVX2(2, ksw_extend_avx2_u8);
|
||||||
#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]);
|
|
||||||
|
|
||||||
#ifdef SHOW_PERF
|
// avx2 heuristics
|
||||||
start_time = get_mseconds();
|
PERFORMANCE_TEST_AVX2(3, ksw_extend_avx2_heuristics);
|
||||||
#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]);
|
|
||||||
|
|
||||||
// 更新query和target位置信息
|
// avx2 u8 heuristics
|
||||||
cur_query_pos += info_arr[i][0];
|
PERFORMANCE_TEST_AVX2(4, ksw_extend_avx2_u8_heuristics);
|
||||||
cur_target_pos += info_arr[i][1];
|
|
||||||
// fprintf(stderr, "%d %d %d %d %d %d %d\n", score_normal, qle, tle, gtle, gscore, max_off[0], max_off[1]);
|
|
||||||
}
|
|
||||||
/*
|
|
||||||
// avx2 sw
|
|
||||||
cur_query_pos = 0;
|
|
||||||
cur_target_pos = 0;
|
|
||||||
for (i = 0; i < block_line_num; ++i)
|
|
||||||
{
|
|
||||||
#ifdef SHOW_PERF
|
|
||||||
int64_t start_time = get_mseconds();
|
|
||||||
#endif
|
|
||||||
score_avx2 += ksw_avx2(
|
|
||||||
info_arr[i][0],
|
|
||||||
(uint8_t *)query_arr + cur_query_pos,
|
|
||||||
info_arr[i][1],
|
|
||||||
(uint8_t *)target_arr + cur_target_pos,
|
|
||||||
0, 5, mat, 6, 1, 6, 1,
|
|
||||||
1, 4,
|
|
||||||
100, 5, 100,
|
|
||||||
info_arr[i][2],
|
|
||||||
&qle, &tle, >le, &gscore, &max_off[0]);
|
|
||||||
#ifdef SHOW_PERF
|
|
||||||
time_sw_avx2 += get_mseconds() - start_time;
|
|
||||||
#endif
|
|
||||||
// 更新query和target位置信息
|
|
||||||
cur_query_pos += info_arr[i][0];
|
|
||||||
cur_target_pos += info_arr[i][1];
|
|
||||||
// fprintf(stderr, "%d %d %d %d %d %d %d\n", score_avx2, qle, tle, gtle, gscore, max_off[0], max_off[1]);
|
|
||||||
}
|
|
||||||
|
|
||||||
// avx2 u8 sw
|
// avx2 mem aligned
|
||||||
cur_query_pos = 0;
|
PERFORMANCE_TEST_AVX2(5, ksw_extend_avx2_aligned);
|
||||||
cur_target_pos = 0;
|
|
||||||
for (i = 0; i < block_line_num; ++i)
|
// avx2 u8 mem aligned
|
||||||
{
|
PERFORMANCE_TEST_AVX2(6, ksw_extend_avx2_u8_aligned);
|
||||||
#ifdef SHOW_PERF
|
|
||||||
int64_t start_time = get_mseconds();
|
|
||||||
#endif
|
|
||||||
score_avx2_u8 += ksw_avx2_u8(
|
|
||||||
info_arr[i][0],
|
|
||||||
(uint8_t *)query_arr + cur_query_pos,
|
|
||||||
info_arr[i][1],
|
|
||||||
(uint8_t *)target_arr + cur_target_pos,
|
|
||||||
0, 5, mat, 6, 1, 6, 1,
|
|
||||||
1, 4,
|
|
||||||
100, 5, 100,
|
|
||||||
info_arr[i][2],
|
|
||||||
&qle, &tle, >le, &gscore, &max_off[0]);
|
|
||||||
#ifdef SHOW_PERF
|
|
||||||
time_sw_avx2_u8 += get_mseconds() - start_time;
|
|
||||||
#endif
|
|
||||||
// 更新query和target位置信息
|
|
||||||
cur_query_pos += info_arr[i][0];
|
|
||||||
cur_target_pos += info_arr[i][1];
|
|
||||||
// fprintf(stderr, "%d %d %d %d %d %d %d\n", score_normal, qle, tle, gtle, gscore, max_off[0], max_off[1]);
|
|
||||||
}
|
|
||||||
*/
|
|
||||||
// fprintf(stderr, "%d %d \n", block_line_num, total_line_num);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// fprintf(stderr, "%d \n", score_normal);
|
|
||||||
|
|
||||||
#ifdef SHOW_PERF
|
#ifdef SHOW_PERF
|
||||||
fprintf(stderr, "time_sw_normal: %f s; score: %d\n", time_sw_normal / DIVIDE_BY, score_normal_total);
|
char *kernel_names[7] = {
|
||||||
fprintf(stderr, "time_bsw_avx2: %f s; score: %d\n", time_bsw_avx2 / DIVIDE_BY, score_bsw_avx2_total);
|
"normal",
|
||||||
fprintf(stderr, "time_sw_avx2: %f s; score: %d\n", time_sw_avx2 / DIVIDE_BY, score_avx2_total);
|
"avx2",
|
||||||
fprintf(stderr, "time_sw_avx2_u8: %f s; score: %d\n", time_sw_avx2_u8 / DIVIDE_BY, score_avx2_u8_total);
|
"avx2_u8",
|
||||||
fprintf(stderr, "time_sw_avx2_u8_aligned: %f s; score: %d\n", time_sw_avx2_u8_aligned / DIVIDE_BY, score_avx2_u8_aligned_total);
|
"avx2_heuristics",
|
||||||
fprintf(stderr, "thread mem capacity: %d\t%d\n", tmem.capacity, tmem_u8.capacity);
|
"avx2_u8_heuristics",
|
||||||
fprintf(stderr, "time_bsw_init: %f s\n", time_bsw_init / DIVIDE_BY);
|
"avx2_aligned",
|
||||||
// fprintf(stderr, "time_bsw_main_loop: %f s\n", (time_bsw_main_loop) / DIVIDE_BY);
|
"avx2_u8_aligned"};
|
||||||
// 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);
|
for (i = 0; i < KERNEL_NUM; ++i)
|
||||||
// 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, "[%18s] time: %9.6f s; score: %d\n", kernel_names[i], time_sw[i] / DIVIDE_BY, score_total[i]);
|
||||||
// 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);
|
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
if (query_f != 0)
|
if (query_f != 0)
|
||||||
|
|
@ -422,6 +263,4 @@ int main(int argc, char *argv[])
|
||||||
fclose(avx2_u8_out_f);
|
fclose(avx2_u8_out_f);
|
||||||
if (normal_out_f != 0)
|
if (normal_out_f != 0)
|
||||||
fclose(normal_out_f);
|
fclose(normal_out_f);
|
||||||
if (bsw_avx2_out_f != 0)
|
|
||||||
fclose(bsw_avx2_out_f);
|
|
||||||
}
|
}
|
||||||
|
|
|
||||||
|
|
@ -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
|
||||||
|
|
@ -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
|
||||||
|
|
@ -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
|
||||||
|
|
@ -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
|
||||||
|
|
@ -9,6 +9,7 @@
|
||||||
|
|
||||||
#include "thread_mem.h"
|
#include "thread_mem.h"
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
|
#include <string.h>
|
||||||
|
|
||||||
// 创建
|
// 创建
|
||||||
thread_mem_t *create_thread_mem()
|
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;
|
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
|
// 将不用的内存归还给thread mem
|
||||||
void thread_mem_release(thread_mem_t *tmem, size_t byte_cnt)
|
void thread_mem_release(thread_mem_t *tmem, size_t byte_cnt)
|
||||||
{
|
{
|
||||||
|
|
|
||||||
|
|
@ -16,7 +16,7 @@
|
||||||
#define MEM_ALIGN_BYTE 8
|
#define MEM_ALIGN_BYTE 8
|
||||||
#define MEM_MOVE_BIT 3
|
#define MEM_MOVE_BIT 3
|
||||||
|
|
||||||
typedef struct
|
typedef struct _thread_mem_t
|
||||||
{
|
{
|
||||||
size_t occupied; // 已经占用的容量（字节数） 对齐的
|
size_t occupied; // 已经占用的容量（字节数） 对齐的
|
||||||
size_t capacity; // 总容量（字节数）
|
size_t capacity; // 总容量（字节数）
|
||||||
|
|
@ -29,9 +29,13 @@ thread_mem_t *create_thread_mem();
|
||||||
void init_thread_mem(thread_mem_t *tmem);
|
void init_thread_mem(thread_mem_t *tmem);
|
||||||
// 初始化并开辟一定量的内存
|
// 初始化并开辟一定量的内存
|
||||||
void thread_mem_init_alloc(thread_mem_t *tmem, size_t byte_cnt);
|
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(thread_mem_t *tmem, size_t byte_cnt);
|
||||||
|
|
||||||
|
// 请求内存并初始化为零
|
||||||
|
void *thread_mem_request_and_clean(thread_mem_t *tmem, size_t byte_cnt);
|
||||||
|
|
||||||
// 将不用的内存归还给thread mem
|
// 将不用的内存归还给thread mem
|
||||||
void thread_mem_release(thread_mem_t *tmem, size_t byte_cnt);
|
void thread_mem_release(thread_mem_t *tmem, size_t byte_cnt);
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -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 <string.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
|
||||||
|
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]];
|
||||||
|
}
|
||||||
|
|
@ -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
|
||||||
Loading…
Reference in New Issue