sc消融实验版本
This commit is contained in:
parent
c04c6541ef
commit
d16f8487ce
|
|
@ -1,3 +1,4 @@
|
|||
*.sam
|
||||
output/
|
||||
input/
|
||||
*.[oa]
|
||||
|
|
|
|||
|
|
@ -11,9 +11,9 @@
|
|||
"request": "launch",
|
||||
"program": "${workspaceRoot}/sw_perf",
|
||||
"args": [
|
||||
"/home/zzh/data/sw/q_s.fa",
|
||||
"/home/zzh/data/sw/t_s.fa",
|
||||
"/home/zzh/data/sw/i_s.txt"
|
||||
"/home/zzh/sw/small/q_s.fa",
|
||||
"/home/zzh/sw/small/t_s.fa",
|
||||
"/home/zzh/sw/small/i_s.txt"
|
||||
],
|
||||
"cwd": "${workspaceFolder}", // 当前工作路径：当前文件所在的工作空间
|
||||
},
|
||||
|
|
|
|||
18
Makefile
18
Makefile
|
|
@ -1,21 +1,17 @@
|
|||
CC= gcc
|
||||
#CFLAGS= -g -Wall -Wno-unused-function -mavx2
|
||||
CFLAGS= -Wall -Wno-unused-function -O2 -mavx2
|
||||
DFLAGS= -DSHOW_PERF -DDEBUG_RETURN_VALUE
|
||||
CFLAGS= -Wall -Wno-unused-function -mavx2 -g -O2
|
||||
DFLAGS= -DSHOW_PERF #-DDEBUG_RETURN_VALUE
|
||||
#DFLAGS= -DSHOW_PERF -DDEBUG_OUT -DDEBUG_RETURN_VALUE
|
||||
PROG= sw_perf
|
||||
PROG2= get_line
|
||||
INCLUDES=
|
||||
LIBS=
|
||||
SUBDIRS= .
|
||||
OBJS= ksw_ext_normal.o \
|
||||
ksw_ext_avx2.o \
|
||||
ksw_ext_avx2_u8.o \
|
||||
ksw_ext_cuda.o \
|
||||
ksw_ext_avx2_heuristics.o \
|
||||
ksw_ext_avx2_u8_heuristics.o \
|
||||
ksw_ext_avx2_aligned.o \
|
||||
ksw_ext_avx2_u8_aligned.o \
|
||||
OBJS= normal.o \
|
||||
normal_pruning.o \
|
||||
avx2_u8.o \
|
||||
avx2_u8_pruning.o \
|
||||
thread_mem.o \
|
||||
utils.o
|
||||
|
||||
|
|
@ -43,4 +39,4 @@ clean:
|
|||
depend:
|
||||
( LC_ALL=C ; export LC_ALL; makedepend -Y -- $(CFLAGS) $(DFLAGS) $(CPPFLAGS) -- *.c )
|
||||
|
||||
# DO NOT DELETE THIS LINE -- make depend depends on it.
|
||||
# DO NOT DELETE THIS LINE -- make depend depends on it.
|
||||
|
|
|
|||
|
|
@ -16,12 +16,20 @@
|
|||
#define UNLIKELY(x) (x)
|
||||
#endif
|
||||
|
||||
#define KSW_EQUAL
|
||||
|
||||
#undef MAX
|
||||
#undef MIN
|
||||
#define MAX(x, y) ((x) > (y) ? (x) : (y))
|
||||
#define MIN(x, y) ((x) < (y) ? (x) : (y))
|
||||
#define SIMD_WIDTH 32
|
||||
|
||||
typedef struct
|
||||
{
|
||||
size_t m;
|
||||
uint8_t *addr;
|
||||
} buf_t;
|
||||
|
||||
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},
|
||||
|
|
@ -58,10 +66,9 @@ static const uint8_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = {
|
|||
|
||||
// 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)
|
||||
// 初始化变量
|
||||
// const int 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; \
|
||||
|
|
@ -77,8 +84,8 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14,
|
|||
oe_ins_vec = _mm256_set1_epi8(oe_ins); \
|
||||
e_del_vec = _mm256_set1_epi8(e_del); \
|
||||
e_ins_vec = _mm256_set1_epi8(e_ins); \
|
||||
__m256i match_sc_vec = _mm256_set1_epi8(base_match_score); \
|
||||
__m256i mis_sc_vec = _mm256_set1_epi8(base_mis_score); \
|
||||
__m256i match_sc_vec = _mm256_set1_epi8(a); \
|
||||
__m256i mis_sc_vec = _mm256_set1_epi8(b); \
|
||||
__m256i amb_sc_vec = _mm256_set1_epi8(1); \
|
||||
__m256i amb_vec = _mm256_set1_epi8(4); \
|
||||
reverse_mask_vec = _mm256_loadu_si256((__m256i *)(reverse_mask)); \
|
||||
|
|
@ -163,8 +170,8 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14,
|
|||
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) \
|
||||
m = MAX(m, maxVal[0]); \
|
||||
if (maxVal[0] > 0) \
|
||||
{ \
|
||||
for (j = beg, i = iend; j <= end; j += SIMD_WIDTH, i -= SIMD_WIDTH) \
|
||||
{ \
|
||||
|
|
@ -196,32 +203,35 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14,
|
|||
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上位置差的 最大值
|
||||
int 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 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上位置差的 最大值)
|
||||
{
|
||||
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 is_left = 0; // 是不是向左扩展
|
||||
int a = 1; // 碱基match时的分数
|
||||
int b = 4; // 碱基mismatch时的惩罚分数（正数）
|
||||
int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH;
|
||||
int i, ibeg, iend, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
int i, ibeg, D, j, k, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
int Dloop = tlen + qlen; // 循环跳出条件
|
||||
int span, beg1, end1; // 边界条件计算
|
||||
int col_size = qlen + 2 + SIMD_WIDTH;
|
||||
|
|
@ -230,27 +240,35 @@ int ksw_extend_avx2_u8(thread_mem_t *tmem,
|
|||
|
||||
SIMD_INIT; // 初始化simd用的数据
|
||||
|
||||
assert(init_score > 0);
|
||||
// fprintf(stdout, "%d\n", h0);
|
||||
|
||||
assert(h0 > 0);
|
||||
|
||||
// allocate memory
|
||||
// mem = malloc(mem_size);
|
||||
mem = thread_mem_request(tmem, mem_size);
|
||||
mem = malloc(mem_size);
|
||||
//if (buf->m < mem_size)
|
||||
//{
|
||||
// buf->m = mem_size;
|
||||
// buf->addr = realloc(buf->addr, mem_size);
|
||||
//}
|
||||
//mem = buf->addr;
|
||||
|
||||
qtmem = &mem[0];
|
||||
seq = (uint8_t *)&qtmem[0];
|
||||
ref = (uint8_t *)&qtmem[seq_size];
|
||||
if (extend_left)
|
||||
if (is_left)
|
||||
{
|
||||
for (i = 0; i < qlen; ++i)
|
||||
seq[i] = query[qlen - 1 - i];
|
||||
for (i = 0; i < tlen; ++i)
|
||||
ref[i + SIMD_WIDTH - 1] = target[tlen - 1 - 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 - 1] = target[i];
|
||||
ref[i + SIMD_WIDTH] = target[i];
|
||||
}
|
||||
|
||||
vmem = &ref[ref_size];
|
||||
|
|
@ -274,73 +292,90 @@ int ksw_extend_avx2_u8(thread_mem_t *tmem,
|
|||
fA1 = &fA[0];
|
||||
fA2 = &fA[col_size];
|
||||
|
||||
// adjust $window_size if it is too large
|
||||
// adjust $w 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 = max_ins > 1 ? max_ins : 1;
|
||||
window_size = window_size < max_ins ? window_size : max_ins;
|
||||
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;
|
||||
window_size = window_size < max_del ? window_size : max_del; // TODO: is this necessary?
|
||||
w = w < max_del ? w : max_del; // TODO: is this necessary?
|
||||
if (tlen < qlen)
|
||||
window_size = MIN(tlen - 1, window_size);
|
||||
w = MIN(tlen - 1, w);
|
||||
|
||||
// DP loop
|
||||
max = init_score, max_i = max_j = -1;
|
||||
max = h0, max_i = max_j = -1;
|
||||
max_ie = -1, gscore = -1;
|
||||
|
||||
;
|
||||
max_off = 0;
|
||||
beg = 1;
|
||||
end = qlen;
|
||||
// init init_score
|
||||
hA0[0] = init_score; // 左上角
|
||||
fA1[1] = MAX(0, init_score - (o_ins + e_ins));
|
||||
eA2[0] = init_score;
|
||||
hA1[1] = fA1[1];
|
||||
// init h0
|
||||
hA0[0] = h0; // 左上角
|
||||
|
||||
if (qlen == 0 || tlen == 0)
|
||||
Dloop = 0; // 防止意外情况
|
||||
if (window_size >= qlen)
|
||||
if (w >= qlen)
|
||||
{
|
||||
max_ie = 0;
|
||||
gscore = 0;
|
||||
}
|
||||
|
||||
int m_last = 0;
|
||||
int iend;
|
||||
#ifdef KSW_EQUAL
|
||||
int midx = 1, icheck = 0, checkspecial = 1;
|
||||
int m3 = 0, m2 = 0, m1 = 0;
|
||||
#endif
|
||||
for (D = 1; LIKELY(D < Dloop); ++D)
|
||||
{
|
||||
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
||||
if (D < tlen)
|
||||
beg1 = 1;
|
||||
if (D > tlen)
|
||||
{
|
||||
span = MIN(Dloop - D, w);
|
||||
beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1);
|
||||
}
|
||||
else
|
||||
beg1 = D - tlen + 1;
|
||||
if (D < qlen)
|
||||
end1 = D; // 闭区间
|
||||
else
|
||||
end1 = qlen;
|
||||
// beg1 = MAX(D - window_size, beg1);
|
||||
// end1 = MIN(D + window_size, end1);
|
||||
// beg = MAX(beg1, beg);
|
||||
// end = MIN(end1, end);
|
||||
// if (beg > end)
|
||||
// break;
|
||||
{
|
||||
span = MIN(D - 1, w);
|
||||
beg1 = MAX(1, ((D - w) / 2) + 1);
|
||||
}
|
||||
end1 = MIN(qlen, beg1 + span);
|
||||
|
||||
beg = beg1;
|
||||
end = end1;
|
||||
if (beg < beg1)
|
||||
beg = beg1;
|
||||
if (end > end1)
|
||||
end = end1;
|
||||
if (beg > end)
|
||||
break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug
|
||||
|
||||
iend = D - beg; // ref开始计算的位置，倒序
|
||||
iend = D - (beg - 1); // ref开始计算的位置，倒序
|
||||
span = end - beg;
|
||||
ibeg = iend - span; // 0开始的ref索引位置
|
||||
ibeg = iend - span - 1; // 0开始的ref索引位置
|
||||
|
||||
// 每一轮需要记录的数据
|
||||
int m = 0, mj = -1, mi = -1;
|
||||
max_vec = zero_vec;
|
||||
|
||||
// 处理左边界
|
||||
// 要处理边界
|
||||
// 左边界 处理f (insert)
|
||||
if (ibeg == 0)
|
||||
{
|
||||
hA1[end] = MAX(0, h0 - (o_ins + e_ins * end));
|
||||
m = hA1[end];
|
||||
}
|
||||
// 上边界
|
||||
if (beg == 1)
|
||||
{
|
||||
hA0[0] = eA2[0];
|
||||
mA1[0] = 0;
|
||||
eA1[0] = MAX(0, init_score - (o_del + e_del * (iend + 1)));
|
||||
hA1[0] = MAX(0, h0 - (o_del + e_del * iend));
|
||||
}
|
||||
else if (D & 1)
|
||||
{
|
||||
hA1[beg - 1] = 0;
|
||||
hA2[beg - 1] = 0;
|
||||
}
|
||||
|
||||
for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH)
|
||||
|
|
@ -368,16 +403,37 @@ int ksw_extend_avx2_u8(thread_mem_t *tmem,
|
|||
// 存储结果
|
||||
SIMD_STORE;
|
||||
}
|
||||
// 处理上边界
|
||||
if (ibeg == 0)
|
||||
{
|
||||
fA2[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1)));
|
||||
hA2[end + 1] = fA2[end + 1];
|
||||
mA2[end + 1] = 0;
|
||||
}
|
||||
|
||||
SIMD_FIND_MAX;
|
||||
|
||||
#ifdef KSW_EQUAL
|
||||
if (hA1[0] < 4 && checkspecial)
|
||||
{ // b == 4
|
||||
if (hA1[0] == 3)
|
||||
{
|
||||
icheck = iend + 1;
|
||||
}
|
||||
else if (midx == 2)
|
||||
{
|
||||
m2 = MAX(m2, hA2[midx - 1]);
|
||||
}
|
||||
else
|
||||
{
|
||||
m2 = MAX(m2, hA2[midx - 1]);
|
||||
m1 = MAX(m1, hA2[midx - 2]);
|
||||
}
|
||||
m3 = MAX(m3, hA2[midx]);
|
||||
midx += 1;
|
||||
if (ibeg > icheck)
|
||||
{
|
||||
if (!m1 || !m2 || !m3)
|
||||
break;
|
||||
else
|
||||
checkspecial = 0;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
// 注意最后跳出循环j的值
|
||||
j = end + 1;
|
||||
|
||||
|
|
@ -386,36 +442,59 @@ int ksw_extend_avx2_u8(thread_mem_t *tmem,
|
|||
max_ie = gscore > hA2[qlen] ? max_ie : ibeg;
|
||||
gscore = gscore > hA2[qlen] ? gscore : hA2[qlen];
|
||||
}
|
||||
//if (m == 0 && m_last == 0)
|
||||
// break; // 一定要注意，斜对角遍历和按列遍历的不同点
|
||||
if (m > max)
|
||||
{
|
||||
max = m, max_i = mi, max_j = mj;
|
||||
max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi);
|
||||
}
|
||||
else if (m == max && max_i >= mi && mj > max_j)
|
||||
{
|
||||
max_i = mi, max_j = mj;
|
||||
max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi);
|
||||
}
|
||||
else if (zdrop > 0)
|
||||
{
|
||||
#if 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;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
// 调整计算的边界
|
||||
// for (j = beg; LIKELY(j <= end); ++j)
|
||||
//{
|
||||
// int has_val = hA1[j - 1] | hA2[j];
|
||||
// if (has_val)
|
||||
// break;
|
||||
//}
|
||||
// {
|
||||
// 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;
|
||||
//}
|
||||
// {
|
||||
// int has_val = hA1[j - 1] | hA2[j];
|
||||
// if (has_val)
|
||||
// break;
|
||||
// else
|
||||
// hA0[j - 1] = 0;
|
||||
// }
|
||||
// end = j + 1 <= qlen ? j + 1 : qlen;
|
||||
beg = 1; end = qlen;
|
||||
|
||||
m_last = m;
|
||||
// swap m, h, e, f
|
||||
SWAP_DATA_POINTER;
|
||||
}
|
||||
|
||||
// free(mem);
|
||||
thread_mem_release(tmem, mem_size);
|
||||
free(mem);
|
||||
if (_qle)
|
||||
*_qle = max_j + 1;
|
||||
if (_tle)
|
||||
|
|
@ -16,12 +16,20 @@
|
|||
#define UNLIKELY(x) (x)
|
||||
#endif
|
||||
|
||||
#define KSW_EQUAL
|
||||
|
||||
#undef MAX
|
||||
#undef MIN
|
||||
#define MAX(x, y) ((x) > (y) ? (x) : (y))
|
||||
#define MIN(x, y) ((x) < (y) ? (x) : (y))
|
||||
#define SIMD_WIDTH 32
|
||||
|
||||
typedef struct
|
||||
{
|
||||
size_t m;
|
||||
uint8_t *addr;
|
||||
} buf_t;
|
||||
|
||||
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},
|
||||
|
|
@ -58,15 +66,13 @@ static const uint8_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = {
|
|||
|
||||
// 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};
|
||||
|
||||
#define permute_mask _MM_SHUFFLE(0, 1, 2, 3)
|
||||
// const int permute_mask = _MM_SHUFFLE(0, 1, 2, 3);
|
||||
// #define permute_mask _MM_SHUFFLE(0, 1, 2, 3)
|
||||
#define permute_mask 27
|
||||
// 初始化变量
|
||||
// 初始化变量
|
||||
#define SIMD_INIT \
|
||||
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
||||
__m256i zero_vec; \
|
||||
__m256i max_vec, last_max_vec = _mm256_set1_epi8(init_score); \
|
||||
__m256i max_vec; \
|
||||
__m256i oe_del_vec; \
|
||||
__m256i oe_ins_vec; \
|
||||
__m256i e_del_vec; \
|
||||
|
|
@ -78,8 +84,8 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14,
|
|||
oe_ins_vec = _mm256_set1_epi8(oe_ins); \
|
||||
e_del_vec = _mm256_set1_epi8(e_del); \
|
||||
e_ins_vec = _mm256_set1_epi8(e_ins); \
|
||||
__m256i match_sc_vec = _mm256_set1_epi8(base_match_score); \
|
||||
__m256i mis_sc_vec = _mm256_set1_epi8(base_mis_score); \
|
||||
__m256i match_sc_vec = _mm256_set1_epi8(a); \
|
||||
__m256i mis_sc_vec = _mm256_set1_epi8(b); \
|
||||
__m256i amb_sc_vec = _mm256_set1_epi8(1); \
|
||||
__m256i amb_vec = _mm256_set1_epi8(4); \
|
||||
reverse_mask_vec = _mm256_loadu_si256((__m256i *)(reverse_mask)); \
|
||||
|
|
@ -152,51 +158,33 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14,
|
|||
mn_vec = _mm256_and_si256(mn_vec, h_vec_mask[end - j]); \
|
||||
hn_vec = _mm256_and_si256(hn_vec, h_vec_mask[end - j]);
|
||||
|
||||
// cmp_max = _mm256_xor_si256(last_max_vec, cmp_max);
|
||||
// last_max_vec = _mm256_set1_epi8(m);
|
||||
// 找最大值和位置
|
||||
#define SIMD_FIND_MAX \
|
||||
__m256i cmp_max = _mm256_max_epu8(max_vec, last_max_vec); \
|
||||
cmp_max = _mm256_xor_si256(cmp_max, last_max_vec); \
|
||||
cmp_max = _mm256_cmpeq_epi8(cmp_max, zero_vec); \
|
||||
uint32_t cmp_result = _mm256_movemask_epi8(cmp_max); \
|
||||
if (cmp_result != 4294967295) \
|
||||
{ \
|
||||
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]; \
|
||||
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; \
|
||||
for (; mj + 1 < qlen && mi + 1 < tlen; mj++, mi++) \
|
||||
{ \
|
||||
if (seq[mj + 1] == ref[mi + 1 + SIMD_WIDTH]) \
|
||||
{ \
|
||||
m += base_match_score; \
|
||||
} \
|
||||
else \
|
||||
{ \
|
||||
break; \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
last_max_vec = _mm256_set1_epi8(m); \
|
||||
#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 = MAX(m, maxVal[0]); \
|
||||
if (maxVal[0] > 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; \
|
||||
} \
|
||||
} \
|
||||
}
|
||||
|
||||
// 每轮迭代后，交换数组
|
||||
|
|
@ -215,64 +203,71 @@ static const uint8_t reverse_mask[SIMD_WIDTH] = {7, 6, 5, 4, 3, 2, 1, 0, 15, 14,
|
|||
mA1 = mA2; \
|
||||
mA2 = tmp;
|
||||
|
||||
// uint8_t mem1[102400];
|
||||
|
||||
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 avx2_u8_pruning(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的惩罚系数
|
||||
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上位置差的 最大值)
|
||||
{
|
||||
// return init_score;
|
||||
uint8_t *mA, *hA, *eA, *fA, *mA1, *mA2, *hA0, *hA1, *eA1, *fA1, *hA2, *eA2, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M
|
||||
uint8_t *seq, *ref;
|
||||
uint8_t *mem, *qtmem, *vmem;
|
||||
int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH;
|
||||
int i, ibeg, iend, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
int i, ibeg, D, j, k, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
int Dloop = tlen + qlen; // 循环跳出条件
|
||||
int span, beg1, end1; // 边界条件计算
|
||||
int col_size = qlen + 2 + SIMD_WIDTH;
|
||||
int val_mem_size = (col_size * 9 + 31) >> 5 << 5; // 32字节的整数倍
|
||||
int mem_size = seq_size + ref_size + val_mem_size;
|
||||
|
||||
int is_left = 0; // 是不是向左扩展
|
||||
int a = 1; // 碱基match时的分数
|
||||
int b = 4; // 碱基mismatch时的惩罚分数（正数）
|
||||
|
||||
SIMD_INIT; // 初始化simd用的数据
|
||||
|
||||
assert(init_score > 0);
|
||||
assert(h0 > 0);
|
||||
|
||||
// allocate memory
|
||||
// mem = malloc(mem_size);
|
||||
mem = thread_mem_request(tmem, mem_size);
|
||||
mem = malloc(mem_size);
|
||||
// if (buf->m < mem_size)
|
||||
//{
|
||||
// buf->m = mem_size;
|
||||
// buf->addr = realloc(buf->addr, mem_size);
|
||||
// }
|
||||
// mem = buf->addr;
|
||||
|
||||
qtmem = &mem[0];
|
||||
seq = (uint8_t *)&qtmem[0];
|
||||
ref = (uint8_t *)&qtmem[seq_size];
|
||||
if (extend_left)
|
||||
if (is_left)
|
||||
{
|
||||
for (i = 0; i < qlen; ++i)
|
||||
seq[i] = query[qlen - 1 - i];
|
||||
for (i = 0; i < tlen; ++i)
|
||||
ref[i + SIMD_WIDTH - 1] = target[tlen - 1 - 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 - 1] = target[i];
|
||||
ref[i + SIMD_WIDTH] = target[i];
|
||||
}
|
||||
|
||||
vmem = &ref[ref_size];
|
||||
|
|
@ -296,74 +291,90 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem,
|
|||
fA1 = &fA[0];
|
||||
fA2 = &fA[col_size];
|
||||
|
||||
// adjust $window_size if it is too large
|
||||
// adjust $w 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 = max_ins > 1 ? max_ins : 1;
|
||||
window_size = window_size < max_ins ? window_size : max_ins;
|
||||
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;
|
||||
window_size = window_size < max_del ? window_size : max_del; // TODO: is this necessary?
|
||||
w = w < max_del ? w : max_del; // TODO: is this necessary?
|
||||
if (tlen < qlen)
|
||||
window_size = MIN(tlen - 1, window_size);
|
||||
w = MIN(tlen - 1, w);
|
||||
|
||||
// DP loop
|
||||
max = init_score, max_i = max_j = -1;
|
||||
max = h0, max_i = max_j = -1;
|
||||
max_ie = -1, gscore = -1;
|
||||
|
||||
;
|
||||
max_off = 0;
|
||||
beg = 1;
|
||||
end = qlen;
|
||||
// init init_score
|
||||
hA0[0] = init_score; // 左上角
|
||||
fA1[1] = MAX(0, init_score - (o_ins + e_ins));
|
||||
eA2[0] = init_score;
|
||||
hA1[1] = fA1[1];
|
||||
// init h0
|
||||
hA0[0] = h0; // 左上角
|
||||
|
||||
if (qlen == 0 || tlen == 0)
|
||||
Dloop = 0; // 防止意外情况
|
||||
if (window_size >= qlen)
|
||||
if (w >= qlen)
|
||||
{
|
||||
max_ie = 0;
|
||||
gscore = 0;
|
||||
}
|
||||
|
||||
int m_last = 0;
|
||||
int iend;
|
||||
#ifdef KSW_EQUAL
|
||||
int midx = 1, icheck = 0, checkspecial = 1;
|
||||
int m3 = 0, m2 = 0, m1 = 0;
|
||||
#endif
|
||||
for (D = 1; LIKELY(D < Dloop); ++D)
|
||||
{
|
||||
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
||||
if (D < tlen)
|
||||
beg1 = 1;
|
||||
if (D > tlen)
|
||||
{
|
||||
span = MIN(Dloop - D, w);
|
||||
beg1 = MAX(D - tlen + 1, ((D - w) / 2) + 1);
|
||||
}
|
||||
else
|
||||
beg1 = D - tlen + 1;
|
||||
if (D < qlen)
|
||||
end1 = D; // 闭区间
|
||||
else
|
||||
end1 = qlen;
|
||||
beg1 = MAX(D - window_size, beg1);
|
||||
end1 = MIN(D + window_size, end1);
|
||||
{
|
||||
span = MIN(D - 1, w);
|
||||
beg1 = MAX(1, ((D - w) / 2) + 1);
|
||||
}
|
||||
end1 = MIN(qlen, beg1 + span);
|
||||
|
||||
beg = MAX(beg1, beg);
|
||||
end = MIN(end1, end);
|
||||
if (beg < beg1)
|
||||
beg = beg1;
|
||||
if (end > end1)
|
||||
end = end1;
|
||||
if (beg > end)
|
||||
break;
|
||||
break; // 不用计算了，直接跳出，否则hA2没有被赋值，里边是上一轮hA0的值，会出bug
|
||||
|
||||
// beg = beg1;
|
||||
// end = end1;
|
||||
|
||||
iend = D - beg; // ref开始计算的位置，倒序
|
||||
iend = D - (beg - 1); // ref开始计算的位置，倒序
|
||||
span = end - beg;
|
||||
ibeg = iend - span; // 0开始的ref索引位置
|
||||
ibeg = iend - span - 1; // 0开始的ref索引位置
|
||||
|
||||
// 每一轮需要记录的数据
|
||||
int m = 0, mj = -1, mi = -1;
|
||||
max_vec = zero_vec;
|
||||
|
||||
// 处理左边界
|
||||
// 要处理边界
|
||||
// 左边界 处理f (insert)
|
||||
if (ibeg == 0)
|
||||
{
|
||||
hA1[end] = MAX(0, h0 - (o_ins + e_ins * end));
|
||||
m = hA1[end];
|
||||
}
|
||||
// 上边界
|
||||
if (beg == 1)
|
||||
{
|
||||
hA0[0] = eA2[0];
|
||||
mA1[0] = 0;
|
||||
eA1[0] = MAX(0, init_score - (o_del + e_del * (iend + 1)));
|
||||
hA1[0] = MAX(0, h0 - (o_del + e_del * iend));
|
||||
}
|
||||
else if (D & 1)
|
||||
{
|
||||
hA1[beg - 1] = 0;
|
||||
hA2[beg - 1] = 0;
|
||||
}
|
||||
|
||||
for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH)
|
||||
|
|
@ -391,16 +402,37 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem,
|
|||
// 存储结果
|
||||
SIMD_STORE;
|
||||
}
|
||||
// 处理上边界
|
||||
if (ibeg == 0)
|
||||
{
|
||||
fA2[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1)));
|
||||
hA2[end + 1] = fA2[end + 1];
|
||||
mA2[end + 1] = 0;
|
||||
}
|
||||
|
||||
SIMD_FIND_MAX;
|
||||
|
||||
#ifdef KSW_EQUAL
|
||||
if (hA1[0] < 4 && checkspecial)
|
||||
{ // b == 4
|
||||
if (hA1[0] == 3)
|
||||
{
|
||||
icheck = iend + 1;
|
||||
}
|
||||
else if (midx == 2)
|
||||
{
|
||||
m2 = MAX(m2, hA2[midx - 1]);
|
||||
}
|
||||
else
|
||||
{
|
||||
m2 = MAX(m2, hA2[midx - 1]);
|
||||
m1 = MAX(m1, hA2[midx - 2]);
|
||||
}
|
||||
m3 = MAX(m3, hA2[midx]);
|
||||
midx += 1;
|
||||
if (ibeg > icheck)
|
||||
{
|
||||
if (!m1 || !m2 || !m3)
|
||||
break;
|
||||
else
|
||||
checkspecial = 0;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
// 注意最后跳出循环j的值
|
||||
j = end + 1;
|
||||
|
||||
|
|
@ -409,12 +441,33 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem,
|
|||
max_ie = gscore > hA2[qlen] ? max_ie : ibeg;
|
||||
gscore = gscore > hA2[qlen] ? gscore : hA2[qlen];
|
||||
}
|
||||
|
||||
if (m == 0 && m_last == 0)
|
||||
break; // 一定要注意，斜对角遍历和按列遍历的不同点
|
||||
if (m > max)
|
||||
{
|
||||
max = m, max_i = mi, max_j = mj;
|
||||
max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi);
|
||||
}
|
||||
else if (m == max && max_i >= mi && mj > max_j)
|
||||
{
|
||||
max_i = mi, max_j = mj;
|
||||
max_off = max_off > abs(mj - mi) ? max_off : abs(mj - mi);
|
||||
}
|
||||
else if (zdrop > 0)
|
||||
{
|
||||
#if 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;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
// 调整计算的边界
|
||||
for (j = beg; LIKELY(j <= end); ++j)
|
||||
|
|
@ -429,20 +482,17 @@ int ksw_extend_avx2_u8_heuristics(thread_mem_t *tmem,
|
|||
int has_val = hA1[j - 1] | hA2[j];
|
||||
if (has_val)
|
||||
break;
|
||||
else
|
||||
hA0[j - 1] = 0;
|
||||
}
|
||||
end = j + 1 <= qlen ? j + 1 : qlen;
|
||||
|
||||
m_last = m;
|
||||
// swap m, h, e, f
|
||||
SWAP_DATA_POINTER;
|
||||
#ifdef SHOW_PERF
|
||||
// time_bsw_adjust_bound += get_mseconds() - start_time;
|
||||
// start_time = get_mseconds();
|
||||
// time_compare += get_mseconds() - start_time;
|
||||
#endif
|
||||
}
|
||||
|
||||
// free(mem);
|
||||
thread_mem_release(tmem, mem_size);
|
||||
free(mem);
|
||||
if (_qle)
|
||||
*_qle = max_j + 1;
|
||||
if (_tle)
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
/*********************************************************************************************
|
||||
Description: Declarations of sw extend functions
|
||||
|
||||
Copyright : All right reserved by NCIC.ICT
|
||||
|
||||
Author : Zhang Zhonghai
|
||||
Date : 2023/08/23
|
||||
***********************************************************************************************/
|
||||
#ifndef __BSW_H
|
||||
#define __BSW_H
|
||||
#include <stdint.h>
|
||||
|
||||
typedef struct _thread_mem_t thread_mem_t;
|
||||
|
||||
int 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);
|
||||
int normal_pruning(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);
|
||||
int avx2_u8(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);
|
||||
int avx2_u8_pruning(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);
|
||||
|
||||
#endif
|
||||
2
common.h
2
common.h
|
|
@ -15,7 +15,7 @@
|
|||
#define MAX(x, y) ((x) > (y) ? (x) : (y))
|
||||
#define MIN(x, y) ((x) < (y) ? (x) : (y))
|
||||
|
||||
#define KERNEL_NUM 7
|
||||
#define KERNEL_NUM 4
|
||||
|
||||
extern FILE *ins_ext_f_arr[KERNEL_NUM],
|
||||
*del_ext_f_arr[KERNEL_NUM],
|
||||
|
|
|
|||
164
ksw_ext.h
164
ksw_ext.h
|
|
@ -1,164 +0,0 @@
|
|||
/*********************************************************************************************
|
||||
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
|
||||
473
ksw_ext_avx2.c
473
ksw_ext_avx2.c
|
|
@ -1,473 +0,0 @@
|
|||
#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"
|
||||
#include "common.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
|
||||
|
||||
#define SIMD_WIDTH 16
|
||||
|
||||
static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = {
|
||||
{0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}};
|
||||
|
||||
// #define permute_mask _MM_SHUFFLE(0, 1, 2, 3)
|
||||
#define permute_mask 27
|
||||
// 初始化变量
|
||||
#define SIMD_INIT \
|
||||
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
||||
__m256i zero_vec; \
|
||||
__m256i max_vec; \
|
||||
__m256i oe_del_vec; \
|
||||
__m256i oe_ins_vec; \
|
||||
__m256i e_del_vec; \
|
||||
__m256i e_ins_vec; \
|
||||
__m256i h_vec_mask[SIMD_WIDTH]; \
|
||||
zero_vec = _mm256_setzero_si256(); \
|
||||
oe_del_vec = _mm256_set1_epi16(-oe_del); \
|
||||
oe_ins_vec = _mm256_set1_epi16(-oe_ins); \
|
||||
e_del_vec = _mm256_set1_epi16(-e_del); \
|
||||
e_ins_vec = _mm256_set1_epi16(-e_ins); \
|
||||
__m256i match_sc_vec = _mm256_set1_epi16(base_match_score); \
|
||||
__m256i mis_sc_vec = _mm256_set1_epi16(-base_mis_score); \
|
||||
__m256i amb_sc_vec = _mm256_set1_epi16(-1); \
|
||||
__m256i amb_vec = _mm256_set1_epi16(4); \
|
||||
for (i = 0; i < SIMD_WIDTH; ++i) \
|
||||
h_vec_mask[i] = _mm256_loadu_si256((__m256i *)(&h_vec_int_mask[i]));
|
||||
|
||||
/*
|
||||
* e 表示当前ref的碱基被删除
|
||||
* f 表示当前seq的碱基插入
|
||||
* m 表示当前碱基匹配（可以相等，也可以不想等）
|
||||
* h 表示最大值
|
||||
*/
|
||||
// load向量化数据
|
||||
#define SIMD_LOAD \
|
||||
__m256i m1 = _mm256_loadu_si256((__m256i *)(&mA1[j])); \
|
||||
__m256i e1 = _mm256_loadu_si256((__m256i *)(&eA1[j])); \
|
||||
__m256i m1j1 = _mm256_loadu_si256((__m256i *)(&mA1[j - 1])); \
|
||||
__m256i f1j1 = _mm256_loadu_si256((__m256i *)(&fA1[j - 1])); \
|
||||
__m256i h0j1 = _mm256_loadu_si256((__m256i *)(&hA0[j - 1])); \
|
||||
__m256i qs_vec = _mm256_loadu_si256((__m256i *)(&seq[j - 1])); \
|
||||
__m256i ts_vec = _mm256_loadu_si256((__m256i *)(&ref[i]));
|
||||
|
||||
// 比对ref和seq的序列，计算罚分
|
||||
#define SIMD_CMP_SEQ \
|
||||
ts_vec = _mm256_permute4x64_epi64(ts_vec, permute_mask); \
|
||||
ts_vec = _mm256_shufflelo_epi16(ts_vec, permute_mask); \
|
||||
ts_vec = _mm256_shufflehi_epi16(ts_vec, permute_mask); \
|
||||
__m256i match_mask_vec = _mm256_cmpeq_epi16(qs_vec, ts_vec); \
|
||||
__m256i mis_score_vec = _mm256_andnot_si256(match_mask_vec, mis_sc_vec); \
|
||||
__m256i score_vec = _mm256_and_si256(match_sc_vec, match_mask_vec); \
|
||||
score_vec = _mm256_or_si256(score_vec, mis_score_vec); \
|
||||
__m256i q_amb_mask_vec = _mm256_cmpeq_epi16(qs_vec, amb_vec); \
|
||||
__m256i t_amb_mask_vec = _mm256_cmpeq_epi16(ts_vec, amb_vec); \
|
||||
__m256i amb_mask_vec = _mm256_or_si256(q_amb_mask_vec, t_amb_mask_vec); \
|
||||
score_vec = _mm256_andnot_si256(amb_mask_vec, score_vec); \
|
||||
__m256i amb_score_vec = _mm256_and_si256(amb_mask_vec, amb_sc_vec); \
|
||||
score_vec = _mm256_or_si256(score_vec, amb_score_vec);
|
||||
|
||||
// 向量化计算h, e, f, m
|
||||
#define SIMD_COMPUTE \
|
||||
__m256i en_vec0 = _mm256_add_epi16(m1, oe_del_vec); \
|
||||
__m256i en_vec1 = _mm256_add_epi16(e1, e_del_vec); \
|
||||
__m256i en_vec = _mm256_max_epi16(en_vec0, en_vec1); \
|
||||
__m256i fn_vec0 = _mm256_add_epi16(m1j1, oe_ins_vec); \
|
||||
__m256i fn_vec1 = _mm256_add_epi16(f1j1, e_ins_vec); \
|
||||
__m256i fn_vec = _mm256_max_epi16(fn_vec0, fn_vec1); \
|
||||
__m256i mn_vec0 = _mm256_add_epi16(h0j1, score_vec); \
|
||||
__m256i mn_mask = _mm256_cmpgt_epi16(h0j1, zero_vec); \
|
||||
__m256i mn_vec = _mm256_and_si256(mn_vec0, mn_mask); \
|
||||
__m256i hn_vec0 = _mm256_max_epi16(en_vec, fn_vec); \
|
||||
__m256i hn_vec = _mm256_max_epi16(hn_vec0, mn_vec); \
|
||||
en_vec = _mm256_max_epi16(en_vec, zero_vec); \
|
||||
fn_vec = _mm256_max_epi16(fn_vec, zero_vec); \
|
||||
mn_vec = _mm256_max_epi16(mn_vec, zero_vec); \
|
||||
hn_vec = _mm256_max_epi16(hn_vec, zero_vec);
|
||||
|
||||
/* int16_t *t_ptr = (int16_t *)&ts_vec; \
|
||||
fprintf(stderr, "D: %d, ibeg: %d, iend: %d, jbeg: %d, jend: %d, %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d \n", \
|
||||
D, ibeg, iend, beg, end, \
|
||||
t_ptr[0], t_ptr[1], t_ptr[2], t_ptr[3], \
|
||||
t_ptr[4], t_ptr[5], t_ptr[6], t_ptr[7], \
|
||||
t_ptr[8], t_ptr[9], t_ptr[10], t_ptr[11], \
|
||||
t_ptr[12], t_ptr[13], t_ptr[14], t_ptr[15]);
|
||||
*/
|
||||
|
||||
// 存储向量化结果
|
||||
#define SIMD_STORE \
|
||||
max_vec = _mm256_max_epi16(max_vec, hn_vec); \
|
||||
_mm256_storeu_si256((__m256i *)&eA2[j], en_vec); \
|
||||
_mm256_storeu_si256((__m256i *)&fA2[j], fn_vec); \
|
||||
_mm256_storeu_si256((__m256i *)&mA2[j], mn_vec); \
|
||||
_mm256_storeu_si256((__m256i *)&hA2[j], hn_vec);
|
||||
|
||||
// 去除多余的部分
|
||||
#define SIMD_REMOVE_EXTRA \
|
||||
en_vec = _mm256_and_si256(en_vec, h_vec_mask[end - j]); \
|
||||
fn_vec = _mm256_and_si256(fn_vec, h_vec_mask[end - j]); \
|
||||
mn_vec = _mm256_and_si256(mn_vec, h_vec_mask[end - j]); \
|
||||
hn_vec = _mm256_and_si256(hn_vec, h_vec_mask[end - j]);
|
||||
|
||||
// 找最大值和位置
|
||||
#define SIMD_FIND_MAX \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 2)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 4)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 6)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 8)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_permute2x128_si256(max_vec, max_vec, 0x01)); \
|
||||
int16_t *maxVal = (int16_t *)&max_vec; \
|
||||
m = maxVal[0]; \
|
||||
if (m > 0) \
|
||||
{ \
|
||||
for (j = beg, i = iend; j <= end; j += SIMD_WIDTH, i -= SIMD_WIDTH) \
|
||||
{ \
|
||||
__m256i h2_vec = _mm256_loadu_si256((__m256i *)(&hA2[j])); \
|
||||
__m256i vcmp = _mm256_cmpeq_epi16(h2_vec, max_vec); \
|
||||
uint32_t mask = _mm256_movemask_epi8(vcmp); \
|
||||
if (mask > 0) \
|
||||
{ \
|
||||
int pos = SIMD_WIDTH - 1 - ((__builtin_clz(mask)) >> 1); \
|
||||
mj = j - 1 + pos; \
|
||||
mi = i - 1 - pos; \
|
||||
} \
|
||||
} \
|
||||
}
|
||||
|
||||
// 每轮迭代后，交换数组
|
||||
#define SWAP_DATA_POINTER \
|
||||
int16_t *tmp = hA0; \
|
||||
hA0 = hA1; \
|
||||
hA1 = hA2; \
|
||||
hA2 = tmp; \
|
||||
tmp = eA1; \
|
||||
eA1 = eA2; \
|
||||
eA2 = tmp; \
|
||||
tmp = fA1; \
|
||||
fA1 = fA2; \
|
||||
fA2 = tmp; \
|
||||
tmp = mA1; \
|
||||
mA1 = mA2; \
|
||||
mA2 = tmp;
|
||||
|
||||
int ksw_extend_avx2(thread_mem_t *tmem,
|
||||
int qlen, // query length 待匹配段碱基的query长度
|
||||
const uint8_t *query, // read碱基序列
|
||||
int tlen, // target length reference的长度
|
||||
const uint8_t *target, // reference序列
|
||||
int extend_left, // 是不是向左扩展
|
||||
int o_del, // deletion 错配开始的惩罚系数
|
||||
int e_del, // deletion extension的惩罚系数
|
||||
int o_ins, // insertion 错配开始的惩罚系数
|
||||
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||
int base_match_score, // 碱基match时的分数
|
||||
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||
int *_gtle, // query全部匹配上的target的长度
|
||||
int *_gscore, // query的端到端匹配得分
|
||||
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
||||
{
|
||||
int16_t *mA, *hA, *eA, *fA, *mA1, *mA2, *hA0, *hA1, *eA1, *fA1, *hA2, *eA2, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M
|
||||
int16_t *seq, *ref;
|
||||
uint8_t *mem;
|
||||
int16_t *qtmem, *vmem;
|
||||
int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH;
|
||||
int i, ibeg, iend, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
int Dloop = tlen + qlen; // 循环跳出条件
|
||||
int span, beg1, end1; // 边界条件计算
|
||||
int col_size = qlen + 2 + SIMD_WIDTH;
|
||||
int val_mem_size = (col_size * 9 * 2 + 31) >> 5 << 5; // 32字节的整数倍
|
||||
int mem_size = (seq_size + ref_size) * 2 + val_mem_size;
|
||||
|
||||
SIMD_INIT; // 初始化simd用的数据
|
||||
|
||||
assert(init_score > 0);
|
||||
|
||||
// allocate memory
|
||||
mem = thread_mem_request(tmem, mem_size);
|
||||
qtmem = (int16_t *)&mem[0];
|
||||
seq = &qtmem[0];
|
||||
ref = &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 - 1] = target[tlen - 1 - i];
|
||||
}
|
||||
else
|
||||
{
|
||||
for (i = 0; i < qlen; ++i)
|
||||
seq[i] = query[i];
|
||||
for (i = 0; i < tlen; ++i)
|
||||
ref[i + SIMD_WIDTH - 1] = target[i];
|
||||
}
|
||||
|
||||
vmem = &ref[ref_size];
|
||||
for (i = 0; i < (val_mem_size >> 1); i += SIMD_WIDTH)
|
||||
{
|
||||
_mm256_storeu_si256((__m256i *)&vmem[i], zero_vec);
|
||||
}
|
||||
hA = &vmem[0];
|
||||
mA = &vmem[col_size * 3];
|
||||
eA = &vmem[col_size * 5];
|
||||
fA = &vmem[col_size * 7];
|
||||
|
||||
hA0 = &hA[0];
|
||||
hA1 = &hA[col_size];
|
||||
hA2 = &hA1[col_size];
|
||||
mA1 = &mA[0];
|
||||
mA2 = &mA[col_size];
|
||||
eA1 = &eA[0];
|
||||
eA2 = &eA[col_size];
|
||||
fA1 = &fA[0];
|
||||
fA2 = &fA[col_size];
|
||||
|
||||
// adjust $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; // 左上角
|
||||
fA1[1] = MAX(0, init_score - (o_ins + e_ins));
|
||||
eA2[0] = init_score;
|
||||
hA1[1] = fA1[1];
|
||||
if (qlen == 0 || tlen == 0)
|
||||
Dloop = 0; // 防止意外情况
|
||||
if (window_size >= qlen)
|
||||
{
|
||||
max_ie = 0;
|
||||
gscore = 0;
|
||||
}
|
||||
// fprintf(stderr, "qlen:%d, tlen:%d\n", qlen, tlen);
|
||||
#ifdef DEBUG_OUT
|
||||
int dii, djj;
|
||||
int16_t ins[tlen + 1][qlen + 2];
|
||||
int16_t del[tlen + 1][qlen + 2];
|
||||
int16_t score[tlen + 1][qlen + 2];
|
||||
for (dii = 0; dii <= tlen; ++dii)
|
||||
{
|
||||
for (djj = 0; djj <= qlen; ++djj)
|
||||
{
|
||||
ins[dii][djj] = del[dii][djj] = score[dii][djj] = 0;
|
||||
}
|
||||
}
|
||||
ins[0][0] = del[0][0] = score[0][0] = init_score;
|
||||
ins[0][1] = MAX(0, init_score - (o_ins + e_ins));
|
||||
del[1][0] = MAX(0, init_score - (o_del + e_del));
|
||||
score[0][1] = ins[0][1];
|
||||
score[1][0] = del[1][0];
|
||||
// fprintf(stderr, "%d %d\n", del[1][0], score[1][0]);
|
||||
#endif
|
||||
|
||||
for (D = 1; LIKELY(D < Dloop); ++D)
|
||||
{
|
||||
if (D < tlen)
|
||||
beg1 = 1;
|
||||
else
|
||||
beg1 = D - tlen + 1;
|
||||
if (D < qlen)
|
||||
end1 = D; // 闭区间
|
||||
else
|
||||
end1 = qlen;
|
||||
// beg1 = MAX(D - window_size, beg1);
|
||||
// end1 = MIN(D + window_size, end1);
|
||||
// beg = MAX(beg1, beg);
|
||||
// end = MIN(end1, end);
|
||||
// if (beg > end)
|
||||
// break;
|
||||
|
||||
beg = beg1;
|
||||
end = end1;
|
||||
|
||||
iend = D - beg; // ref开始计算的位置，倒序
|
||||
span = end - beg;
|
||||
ibeg = iend - span; // 0开始的ref索引位置
|
||||
|
||||
// fprintf(stderr, "D:%d, jbeg:%d, jend:%d, ibeg:%d, iend:%d\n", D, beg, end, ibeg, iend);
|
||||
|
||||
// 每一轮需要记录的数据
|
||||
int m = 0, mj = -1, mi = -1;
|
||||
max_vec = zero_vec;
|
||||
|
||||
// 处理左边界
|
||||
if (beg == 1)
|
||||
{
|
||||
hA0[0] = eA2[0];
|
||||
mA1[0] = 0;
|
||||
eA1[0] = MAX(0, init_score - (o_del + e_del * (iend + 1)));
|
||||
#ifdef DEBUG_OUT
|
||||
del[iend + 1][0] = eA1[0];
|
||||
score[iend + 1][0] = eA1[0];
|
||||
#endif
|
||||
}
|
||||
#ifdef DEBUG_OUT
|
||||
// fprintf(stderr, "eA1: %d\n", eA1[0]);
|
||||
// for (djj = beg - 1; djj < end; ++djj)
|
||||
//{
|
||||
// fprintf(stderr, "%d ", hA0[djj]);
|
||||
//}
|
||||
// fprintf(stderr, "\n");
|
||||
#endif
|
||||
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;
|
||||
}
|
||||
// 处理上边界
|
||||
if (ibeg == 0)
|
||||
{
|
||||
fA2[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1)));
|
||||
hA2[end + 1] = fA2[end + 1];
|
||||
mA2[end + 1] = 0;
|
||||
#ifdef DEBUG_OUT
|
||||
ins[0][end + 1] = fA2[end + 1];
|
||||
score[0][end + 1] = fA2[end + 1];
|
||||
#endif
|
||||
}
|
||||
|
||||
SIMD_FIND_MAX;
|
||||
|
||||
#ifdef DEBUG_OUT
|
||||
for (djj = beg; djj <= end; ++djj)
|
||||
{
|
||||
dii = D - djj + 1;
|
||||
// fprintf(stderr, "dii:%d, djj:%d, ", dii, djj);
|
||||
ins[dii][djj] = fA2[djj];
|
||||
del[dii][djj] = eA2[djj];
|
||||
score[dii][djj] = hA2[djj];
|
||||
}
|
||||
// fprintf(stderr, "\n");
|
||||
// fprintf(stderr, "%d, %d\n", hA2[0], hA2[1]);
|
||||
#endif
|
||||
|
||||
// 注意最后跳出循环j的值
|
||||
j = end + 1;
|
||||
|
||||
if (j == qlen + 1)
|
||||
{
|
||||
max_ie = gscore > hA2[qlen] ? max_ie : ibeg;
|
||||
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;
|
||||
//}
|
||||
// end = j + 1 <= qlen ? j + 1 : qlen;
|
||||
// beg = 0;
|
||||
// end = qlen; // uncomment this line for debugging
|
||||
|
||||
// swap m, h, e, f
|
||||
SWAP_DATA_POINTER;
|
||||
}
|
||||
|
||||
#ifdef DEBUG_OUT
|
||||
for (dii = 0; dii <= tlen; ++dii)
|
||||
{
|
||||
for (djj = 0; djj <= qlen; ++djj)
|
||||
{
|
||||
fprintf(score_f_arr[1], "%-4d", score[dii][djj]);
|
||||
fprintf(ins_ext_f_arr[1], "%-4d", ins[dii][djj]);
|
||||
fprintf(del_ext_f_arr[1], "%-4d", del[dii][djj]);
|
||||
}
|
||||
fprintf(score_f_arr[1], "\n");
|
||||
fprintf(ins_ext_f_arr[1], "\n");
|
||||
fprintf(del_ext_f_arr[1], "\n");
|
||||
}
|
||||
#endif
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
|
@ -1,445 +0,0 @@
|
|||
#include <stdlib.h>
|
||||
#include <stdint.h>
|
||||
#include <assert.h>
|
||||
#include <stdio.h>
|
||||
#include <immintrin.h>
|
||||
#include <emmintrin.h>
|
||||
#include "thread_mem.h"
|
||||
#include "common.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 16
|
||||
#define BASE_BYTES 2
|
||||
#define SCORE_BYTES 2
|
||||
#define BOUNDARY_SCORE_NUM 2
|
||||
#define TMP_SCORE_ARRAY_NUM 9
|
||||
#define MEM_ALIGN_BYTES 32
|
||||
#define ALIGN_SHIFT_BITS 5
|
||||
#define SIMD_BYTES 32
|
||||
|
||||
#define AMBIGUOUS_BASE_CODE 4
|
||||
#define AMBIGUOUS_BASE_SCORE -1
|
||||
|
||||
// 32字节对齐（256位）
|
||||
#define align_mem(x) (((x) + 31) >> 5 << 5)
|
||||
#define align_number(x) align_mem(x)
|
||||
|
||||
/* 去掉多余计算的值 */
|
||||
static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = {
|
||||
{0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}};
|
||||
|
||||
// #define permute_mask _MM_SHUFFLE(0, 1, 2, 3)
|
||||
#define permute_mask 27
|
||||
// 初始化变量
|
||||
#define SIMD_INIT \
|
||||
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
||||
__m256i zero_vec; \
|
||||
__m256i max_vec, last_max_vec = _mm256_set1_epi16(init_score); \
|
||||
__m256i oe_del_vec; \
|
||||
__m256i oe_ins_vec; \
|
||||
__m256i e_del_vec; \
|
||||
__m256i e_ins_vec; \
|
||||
__m256i h_vec_mask[SIMD_WIDTH]; \
|
||||
zero_vec = _mm256_setzero_si256(); \
|
||||
oe_del_vec = _mm256_set1_epi16(-oe_del); \
|
||||
oe_ins_vec = _mm256_set1_epi16(-oe_ins); \
|
||||
e_del_vec = _mm256_set1_epi16(-e_del); \
|
||||
e_ins_vec = _mm256_set1_epi16(-e_ins); \
|
||||
__m256i match_sc_vec = _mm256_set1_epi16(base_match_score); \
|
||||
__m256i mis_sc_vec = _mm256_set1_epi16(-base_mis_score); \
|
||||
__m256i amb_sc_vec = _mm256_set1_epi16(AMBIGUOUS_BASE_SCORE); \
|
||||
__m256i amb_vec = _mm256_set1_epi16(AMBIGUOUS_BASE_CODE); \
|
||||
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])); \
|
||||
__m256i ts_vec = _mm256_loadu_si256((__m256i *)(&ref[i]));
|
||||
|
||||
// 比对ref和seq的序列，计算罚分
|
||||
#define SIMD_CMP_SEQ \
|
||||
/* 将待比对的target序列逆序排列 */ \
|
||||
ts_vec = _mm256_permute4x64_epi64(ts_vec, permute_mask); \
|
||||
ts_vec = _mm256_shufflelo_epi16(ts_vec, permute_mask); \
|
||||
ts_vec = _mm256_shufflehi_epi16(ts_vec, permute_mask); \
|
||||
__m256i match_mask_vec = _mm256_cmpeq_epi16(qs_vec, ts_vec); /* 比对query和target字符序列 */ \
|
||||
__m256i mis_score_vec = _mm256_andnot_si256(match_mask_vec, mis_sc_vec); /* 未匹配上的位置赋值mismatch分数 */ \
|
||||
__m256i score_vec = _mm256_and_si256(match_sc_vec, match_mask_vec); /* 匹配上的位置赋值match分数 */ \
|
||||
score_vec = _mm256_or_si256(score_vec, mis_score_vec); \
|
||||
/* 计算模棱两可的字符（N)的位置的分数 */ \
|
||||
__m256i q_amb_mask_vec = _mm256_cmpeq_epi16(qs_vec, amb_vec); \
|
||||
__m256i t_amb_mask_vec = _mm256_cmpeq_epi16(ts_vec, amb_vec); \
|
||||
__m256i amb_mask_vec = _mm256_or_si256(q_amb_mask_vec, t_amb_mask_vec); \
|
||||
score_vec = _mm256_andnot_si256(amb_mask_vec, score_vec); \
|
||||
__m256i amb_score_vec = _mm256_and_si256(amb_mask_vec, amb_sc_vec); \
|
||||
score_vec = _mm256_or_si256(score_vec, amb_score_vec);
|
||||
|
||||
// 向量化计算h, e, f, m
|
||||
#define SIMD_COMPUTE \
|
||||
__m256i en_vec0 = _mm256_add_epi16(m1, oe_del_vec); \
|
||||
__m256i en_vec1 = _mm256_add_epi16(e1, e_del_vec); \
|
||||
__m256i en_vec = _mm256_max_epi16(en_vec0, en_vec1); \
|
||||
__m256i fn_vec0 = _mm256_add_epi16(m1j1, oe_ins_vec); \
|
||||
__m256i fn_vec1 = _mm256_add_epi16(f1j1, e_ins_vec); \
|
||||
__m256i fn_vec = _mm256_max_epi16(fn_vec0, fn_vec1); \
|
||||
__m256i mn_vec0 = _mm256_add_epi16(h0j1, score_vec); \
|
||||
__m256i mn_mask = _mm256_cmpgt_epi16(h0j1, zero_vec); \
|
||||
__m256i mn_vec = _mm256_and_si256(mn_vec0, mn_mask); \
|
||||
__m256i hn_vec0 = _mm256_max_epi16(en_vec, fn_vec); \
|
||||
__m256i hn_vec = _mm256_max_epi16(hn_vec0, mn_vec); \
|
||||
en_vec = _mm256_max_epi16(en_vec, zero_vec); \
|
||||
fn_vec = _mm256_max_epi16(fn_vec, zero_vec); \
|
||||
mn_vec = _mm256_max_epi16(mn_vec, zero_vec); \
|
||||
hn_vec = _mm256_max_epi16(hn_vec, zero_vec);
|
||||
|
||||
#define SIMD_STORE \
|
||||
max_vec = _mm256_max_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 \
|
||||
__m256i cmp_max = _mm256_cmpgt_epi16(max_vec, last_max_vec); \
|
||||
uint32_t cmp_result = _mm256_movemask_epi8(cmp_max); \
|
||||
if (cmp_result > 0) \
|
||||
{ \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 2)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 4)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 6)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 8)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_permute2x128_si256(max_vec, max_vec, 0x01)); \
|
||||
int16_t *maxVal = (int16_t *)&max_vec; \
|
||||
m = maxVal[0]; \
|
||||
for (j = aligned_read_start_pos, i = aligned_ref_end_pos; j <= end; j += SIMD_WIDTH, i -= SIMD_WIDTH) \
|
||||
{ \
|
||||
__m256i h2_vec = _mm256_loadu_si256((__m256i *)(&hA2[j])); \
|
||||
__m256i vcmp = _mm256_cmpeq_epi16(h2_vec, max_vec); \
|
||||
uint32_t mask = _mm256_movemask_epi8(vcmp); \
|
||||
if (mask > 0) \
|
||||
{ \
|
||||
int pos = SIMD_WIDTH - 1 - ((__builtin_clz(mask)) >> 1); \
|
||||
mj = j - 1 + pos; \
|
||||
mi = i - 1 - pos; \
|
||||
for (; mj + 1 < qlen && mi + 1 < tlen; mj++, mi++) \
|
||||
{ \
|
||||
if (seq[mj + 2] == ref[mi + 1 + SIMD_WIDTH]) \
|
||||
{ \
|
||||
m += base_match_score; \
|
||||
} \
|
||||
else \
|
||||
{ \
|
||||
break; \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
last_max_vec = _mm256_set1_epi16(m); \
|
||||
}
|
||||
|
||||
// 每轮迭代后，交换数组
|
||||
#define SWAP_DATA_POINTER \
|
||||
int16_t *tmp = hA0; \
|
||||
hA0 = hA1; \
|
||||
hA1 = hA2; \
|
||||
hA2 = tmp; \
|
||||
tmp = eA1; \
|
||||
eA1 = eA2; \
|
||||
eA2 = tmp; \
|
||||
tmp = fA1; \
|
||||
fA1 = fA2; \
|
||||
fA2 = tmp; \
|
||||
tmp = mA1; \
|
||||
mA1 = mA2; \
|
||||
mA2 = tmp;
|
||||
|
||||
int ksw_extend_avx2_aligned(thread_mem_t *tmem,
|
||||
int qlen, // query length 待匹配段碱基的query长度
|
||||
const uint8_t *query, // read碱基序列
|
||||
int tlen, // target length reference的长度
|
||||
const uint8_t *target, // reference序列
|
||||
int extend_left, // 是不是向左扩展
|
||||
int o_del, // deletion 错配开始的惩罚系数
|
||||
int e_del, // deletion extension的惩罚系数
|
||||
int o_ins, // insertion 错配开始的惩罚系数
|
||||
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||
int base_match_score, // 碱基match时的分数
|
||||
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||
int *_gtle, // query全部匹配上的target的长度
|
||||
int *_gscore, // query的端到端匹配得分
|
||||
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
||||
{
|
||||
int16_t *mA1, *mA2,
|
||||
*hA0, *hA1, *hA2,
|
||||
*eA1, *eA2,
|
||||
*fA1, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M
|
||||
int16_t *seq, *ref;
|
||||
uint8_t *mem_addr;
|
||||
|
||||
int read_size = align_number(qlen * BASE_BYTES + MEM_ALIGN_BYTES);
|
||||
int ref_size = align_number((tlen + SIMD_WIDTH) * BASE_BYTES);
|
||||
int back_diagnal_num = tlen + qlen; // 循环跳出条件 D从1开始遍历
|
||||
int score_array_size = align_number((qlen + BOUNDARY_SCORE_NUM) * SCORE_BYTES);
|
||||
int score_element_num = score_array_size / SCORE_BYTES;
|
||||
int score_mem_size = score_array_size * TMP_SCORE_ARRAY_NUM;
|
||||
int request_mem_size = read_size + ref_size + score_mem_size + MEM_ALIGN_BYTES * 3; // 左侧内存地址对齐 + 数据向左偏移一个元素 + 末尾SIMD补齐
|
||||
|
||||
int i, ibeg, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
int span, beg1, end1; // 边界条件计算
|
||||
int aligned_read_start_pos, aligned_ref_end_pos;
|
||||
int iend;
|
||||
|
||||
SIMD_INIT; // 初始化simd用的数据
|
||||
|
||||
assert(init_score > 0);
|
||||
|
||||
// allocate memory
|
||||
mem_addr = thread_mem_request(tmem, request_mem_size);
|
||||
mem_addr = (void *)align_mem((uint64_t)mem_addr);
|
||||
ref = (int16_t *)&mem_addr[0];
|
||||
seq = (int16_t *)(mem_addr + ref_size + SIMD_BYTES - BASE_BYTES);
|
||||
if (extend_left)
|
||||
{
|
||||
for (i = 0; i < qlen; ++i)
|
||||
seq[i + 1] = query[qlen - 1 - i];
|
||||
for (i = 0; i < tlen; ++i)
|
||||
ref[i + SIMD_WIDTH - 1] = target[tlen - 1 - i];
|
||||
}
|
||||
else
|
||||
{
|
||||
for (i = 0; i < qlen; ++i)
|
||||
seq[i + 1] = query[i];
|
||||
for (i = 0; i < tlen; ++i)
|
||||
ref[i + SIMD_WIDTH - 1] = target[i];
|
||||
}
|
||||
|
||||
mem_addr += read_size + ref_size + (SIMD_BYTES - SCORE_BYTES);
|
||||
for (i = 0; i < score_mem_size; i += SIMD_BYTES)
|
||||
{
|
||||
_mm256_storeu_si256((__m256i *)&mem_addr[i], zero_vec);
|
||||
}
|
||||
|
||||
hA0 = (int16_t *)&mem_addr[0];
|
||||
hA1 = &hA0[score_element_num];
|
||||
hA2 = &hA1[score_element_num];
|
||||
mA1 = &hA2[score_element_num];
|
||||
mA2 = &mA1[score_element_num];
|
||||
eA1 = &mA2[score_element_num];
|
||||
eA2 = &eA1[score_element_num];
|
||||
fA1 = &eA2[score_element_num];
|
||||
fA2 = &fA1[score_element_num];
|
||||
|
||||
// adjust $window_size if it is too large
|
||||
// get the max score
|
||||
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; // 左上角
|
||||
fA1[1] = MAX(0, init_score - (o_ins + e_ins));
|
||||
eA2[0] = init_score;
|
||||
hA1[1] = fA1[1];
|
||||
|
||||
if (qlen == 0 || tlen == 0)
|
||||
back_diagnal_num = 0; // 防止意外情况
|
||||
if (window_size >= qlen)
|
||||
{
|
||||
max_ie = 0;
|
||||
gscore = 0;
|
||||
}
|
||||
|
||||
for (D = 1; LIKELY(D < back_diagnal_num); ++D)
|
||||
{
|
||||
if (D < tlen)
|
||||
beg1 = 1;
|
||||
else
|
||||
beg1 = D - tlen + 1;
|
||||
if (D < qlen)
|
||||
end1 = D; // 闭区间
|
||||
else
|
||||
end1 = qlen;
|
||||
beg1 = MAX(D - window_size, beg1);
|
||||
end1 = MIN(D + window_size, end1);
|
||||
|
||||
beg = MAX(beg1, beg);
|
||||
end = MIN(end1, end);
|
||||
if (beg > end)
|
||||
break;
|
||||
|
||||
// beg = beg1;
|
||||
// end = end1;
|
||||
|
||||
iend = D - beg; // ref开始计算的位置，倒序
|
||||
span = end - beg;
|
||||
ibeg = iend - span; // 0开始的ref索引位置
|
||||
|
||||
// 每一轮需要记录的数据
|
||||
int m = 0, mj = -1, mi = -1;
|
||||
max_vec = zero_vec;
|
||||
// 处理左边界
|
||||
if (beg == 1)
|
||||
{
|
||||
hA0[0] = eA2[0];
|
||||
mA1[0] = 0;
|
||||
eA1[0] = MAX(0, init_score - (o_del + e_del * (iend + 1)));
|
||||
}
|
||||
|
||||
// aligned_read_start_pos = (beg >> ALIGN_SHIFT_BITS << ALIGN_SHIFT_BITS) + 1;
|
||||
// aligned_ref_end_pos = iend + (beg - aligned_read_start_pos);
|
||||
|
||||
aligned_read_start_pos = beg;
|
||||
aligned_ref_end_pos = iend;
|
||||
|
||||
// fprintf(stderr, "%d\t%d\n", beg, aligned_read_start_pos);
|
||||
|
||||
for (j = aligned_read_start_pos, i = aligned_ref_end_pos; 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;
|
||||
}
|
||||
// 处理上边界
|
||||
if (ibeg == 0)
|
||||
{
|
||||
fA2[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1)));
|
||||
hA2[end + 1] = fA2[end + 1];
|
||||
mA2[end + 1] = 0;
|
||||
}
|
||||
|
||||
SIMD_FIND_MAX;
|
||||
|
||||
// 注意最后跳出循环j的值
|
||||
j = end + 1;
|
||||
|
||||
if (j == qlen + 1) // 遍历到了query最后一个碱基，此时next_max_arr[qlen]为全局匹配的最大分值
|
||||
{
|
||||
max_ie = gscore > hA2[qlen] ? max_ie : ibeg;
|
||||
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;
|
||||
}
|
||||
}
|
||||
end = j + 1 <= qlen ? j + 1 : qlen;
|
||||
// swap m, h, e, f
|
||||
SWAP_DATA_POINTER;
|
||||
}
|
||||
|
||||
thread_mem_release(tmem, request_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;
|
||||
}
|
||||
|
|
@ -1,431 +0,0 @@
|
|||
#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"
|
||||
#include "common.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 16
|
||||
|
||||
#define AMBIGUOUS_BASE_CODE 4
|
||||
#define AMBIGUOUS_BASE_SCORE -1
|
||||
|
||||
/* 去掉多余计算的值 */
|
||||
static const uint16_t h_vec_int_mask[SIMD_WIDTH][SIMD_WIDTH] = {
|
||||
{0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0},
|
||||
{0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}};
|
||||
|
||||
// #define permute_mask _MM_SHUFFLE(0, 1, 2, 3)
|
||||
#define permute_mask 27
|
||||
// 初始化变量
|
||||
#define SIMD_INIT \
|
||||
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
||||
__m256i zero_vec; \
|
||||
__m256i max_vec, last_max_vec = _mm256_set1_epi16(init_score); \
|
||||
__m256i oe_del_vec; \
|
||||
__m256i oe_ins_vec; \
|
||||
__m256i e_del_vec; \
|
||||
__m256i e_ins_vec; \
|
||||
__m256i h_vec_mask[SIMD_WIDTH]; \
|
||||
zero_vec = _mm256_setzero_si256(); \
|
||||
oe_del_vec = _mm256_set1_epi16(-oe_del); \
|
||||
oe_ins_vec = _mm256_set1_epi16(-oe_ins); \
|
||||
e_del_vec = _mm256_set1_epi16(-e_del); \
|
||||
e_ins_vec = _mm256_set1_epi16(-e_ins); \
|
||||
__m256i match_sc_vec = _mm256_set1_epi16(base_match_score); \
|
||||
__m256i mis_sc_vec = _mm256_set1_epi16(-base_mis_score); \
|
||||
__m256i amb_sc_vec = _mm256_set1_epi16(AMBIGUOUS_BASE_SCORE); \
|
||||
__m256i amb_vec = _mm256_set1_epi16(AMBIGUOUS_BASE_CODE); \
|
||||
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 *)(&cur_match_arr[j])); \
|
||||
__m256i e1 = _mm256_loadu_si256((__m256i *)(&cur_del_arr[j])); \
|
||||
__m256i m1j1 = _mm256_loadu_si256((__m256i *)(&cur_match_arr[j - 1])); \
|
||||
__m256i f1j1 = _mm256_loadu_si256((__m256i *)(&cur_ins_arr[j - 1])); \
|
||||
__m256i h0j1 = _mm256_loadu_si256((__m256i *)(&last_max_arr[j - 1])); \
|
||||
__m256i qs_vec = _mm256_loadu_si256((__m256i *)(&seq[j - 1])); \
|
||||
__m256i ts_vec = _mm256_loadu_si256((__m256i *)(&ref[i]));
|
||||
|
||||
// 比对ref和seq的序列，计算罚分
|
||||
#define SIMD_CMP_SEQ \
|
||||
/* 将待比对的target序列逆序排列 */ \
|
||||
ts_vec = _mm256_permute4x64_epi64(ts_vec, permute_mask); \
|
||||
ts_vec = _mm256_shufflelo_epi16(ts_vec, permute_mask); \
|
||||
ts_vec = _mm256_shufflehi_epi16(ts_vec, permute_mask); \
|
||||
__m256i match_mask_vec = _mm256_cmpeq_epi16(qs_vec, ts_vec); /* 比对query和target字符序列 */ \
|
||||
__m256i mis_score_vec = _mm256_andnot_si256(match_mask_vec, mis_sc_vec); /* 未匹配上的位置赋值mismatch分数 */ \
|
||||
__m256i score_vec = _mm256_and_si256(match_sc_vec, match_mask_vec); /* 匹配上的位置赋值match分数 */ \
|
||||
score_vec = _mm256_or_si256(score_vec, mis_score_vec); \
|
||||
/* 计算模棱两可的字符（N)的位置的分数 */ \
|
||||
__m256i q_amb_mask_vec = _mm256_cmpeq_epi16(qs_vec, amb_vec); \
|
||||
__m256i t_amb_mask_vec = _mm256_cmpeq_epi16(ts_vec, amb_vec); \
|
||||
__m256i amb_mask_vec = _mm256_or_si256(q_amb_mask_vec, t_amb_mask_vec); \
|
||||
score_vec = _mm256_andnot_si256(amb_mask_vec, score_vec); \
|
||||
__m256i amb_score_vec = _mm256_and_si256(amb_mask_vec, amb_sc_vec); \
|
||||
score_vec = _mm256_or_si256(score_vec, amb_score_vec);
|
||||
|
||||
// 向量化计算h, e, f, m
|
||||
#define SIMD_COMPUTE \
|
||||
__m256i en_vec0 = _mm256_add_epi16(m1, oe_del_vec); \
|
||||
__m256i en_vec1 = _mm256_add_epi16(e1, e_del_vec); \
|
||||
__m256i en_vec = _mm256_max_epi16(en_vec0, en_vec1); \
|
||||
__m256i fn_vec0 = _mm256_add_epi16(m1j1, oe_ins_vec); \
|
||||
__m256i fn_vec1 = _mm256_add_epi16(f1j1, e_ins_vec); \
|
||||
__m256i fn_vec = _mm256_max_epi16(fn_vec0, fn_vec1); \
|
||||
__m256i mn_vec0 = _mm256_add_epi16(h0j1, score_vec); \
|
||||
__m256i mn_mask = _mm256_cmpgt_epi16(h0j1, zero_vec); \
|
||||
__m256i mn_vec = _mm256_and_si256(mn_vec0, mn_mask); \
|
||||
__m256i hn_vec0 = _mm256_max_epi16(en_vec, fn_vec); \
|
||||
__m256i hn_vec = _mm256_max_epi16(hn_vec0, mn_vec); \
|
||||
en_vec = _mm256_max_epi16(en_vec, zero_vec); \
|
||||
fn_vec = _mm256_max_epi16(fn_vec, zero_vec); \
|
||||
mn_vec = _mm256_max_epi16(mn_vec, zero_vec); \
|
||||
hn_vec = _mm256_max_epi16(hn_vec, zero_vec);
|
||||
|
||||
// 存储向量化结果
|
||||
#define SIMD_STORE \
|
||||
max_vec = _mm256_max_epu8(max_vec, hn_vec); \
|
||||
_mm256_storeu_si256((__m256i *)&next_del_arr[j], en_vec); \
|
||||
_mm256_storeu_si256((__m256i *)&next_ins_arr[j], fn_vec); \
|
||||
_mm256_storeu_si256((__m256i *)&next_match_arr[j], mn_vec); \
|
||||
_mm256_storeu_si256((__m256i *)&next_max_arr[j], hn_vec);
|
||||
|
||||
// 去除多余的部分
|
||||
#define SIMD_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]);
|
||||
|
||||
// 找最大值和位置
|
||||
// last_max_vec = max_vec;
|
||||
#define SIMD_FIND_MAX \
|
||||
__m256i cmp_max = _mm256_cmpgt_epi16(max_vec, last_max_vec); \
|
||||
uint32_t cmp_result = _mm256_movemask_epi8(cmp_max); \
|
||||
if (cmp_result > 0) \
|
||||
{ \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 2)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 4)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 6)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_alignr_epi8(max_vec, max_vec, 8)); \
|
||||
max_vec = _mm256_max_epu16(max_vec, _mm256_permute2x128_si256(max_vec, max_vec, 0x01)); \
|
||||
int16_t *maxVal = (int16_t *)&max_vec; \
|
||||
m = maxVal[0]; \
|
||||
for (j = beg, i = iend; j <= end; j += SIMD_WIDTH, i -= SIMD_WIDTH) \
|
||||
{ \
|
||||
__m256i h2_vec = _mm256_loadu_si256((__m256i *)(&next_max_arr[j])); \
|
||||
__m256i vcmp = _mm256_cmpeq_epi16(h2_vec, max_vec); \
|
||||
uint32_t mask = _mm256_movemask_epi8(vcmp); \
|
||||
if (mask > 0) \
|
||||
{ \
|
||||
int pos = SIMD_WIDTH - 1 - ((__builtin_clz(mask)) >> 1); \
|
||||
mj = j - 1 + pos; \
|
||||
mi = i - 1 - pos; \
|
||||
for (; mj + 1 < qlen && mi + 1 < tlen; mj++, mi++) \
|
||||
{ \
|
||||
if (seq[mj + 1] == ref[mi + 1 + SIMD_WIDTH]) \
|
||||
{ \
|
||||
m += base_match_score; \
|
||||
} \
|
||||
else \
|
||||
{ \
|
||||
break; \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
last_max_vec = _mm256_set1_epi16(m); \
|
||||
}
|
||||
|
||||
// 每轮迭代后，交换数组
|
||||
#define SWAP_DATA_POINTER \
|
||||
int16_t *tmp = last_max_arr; \
|
||||
last_max_arr = cur_max_arr; \
|
||||
cur_max_arr = next_max_arr; \
|
||||
next_max_arr = tmp; \
|
||||
tmp = cur_del_arr; \
|
||||
cur_del_arr = next_del_arr; \
|
||||
next_del_arr = tmp; \
|
||||
tmp = cur_ins_arr; \
|
||||
cur_ins_arr = next_ins_arr; \
|
||||
next_ins_arr = tmp; \
|
||||
tmp = cur_match_arr; \
|
||||
cur_match_arr = next_match_arr; \
|
||||
next_match_arr = tmp;
|
||||
|
||||
int ksw_extend_avx2_heuristics(thread_mem_t *tmem,
|
||||
int qlen, // query length 待匹配段碱基的query长度
|
||||
const uint8_t *query, // read碱基序列
|
||||
int tlen, // target length reference的长度
|
||||
const uint8_t *target, // reference序列
|
||||
int extend_left, // 是不是向左扩展
|
||||
int o_del, // deletion 错配开始的惩罚系数
|
||||
int e_del, // deletion extension的惩罚系数
|
||||
int o_ins, // insertion 错配开始的惩罚系数
|
||||
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||
int base_match_score, // 碱基match时的分数
|
||||
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||
int *_gtle, // query全部匹配上的target的长度
|
||||
int *_gscore, // query的端到端匹配得分
|
||||
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
||||
{
|
||||
// return 0;
|
||||
int16_t *mA, *eA, *hA, *fA,
|
||||
*cur_match_arr, *next_match_arr,
|
||||
*last_max_arr, *cur_max_arr, *next_max_arr,
|
||||
*cur_del_arr, *next_del_arr,
|
||||
*cur_ins_arr, *next_ins_arr; // hA0保存上上个col的H，其他的保存上个H E F M
|
||||
int16_t *seq, *ref;
|
||||
uint8_t *mem;
|
||||
int16_t *qtmem, *vmem;
|
||||
int seq_size = qlen + SIMD_WIDTH, ref_size = tlen + SIMD_WIDTH;
|
||||
int i, ibeg, iend, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
int Dloop = tlen + qlen; // 循环跳出条件 D从1开始遍历
|
||||
int span, beg1, end1; // 边界条件计算
|
||||
int col_size = qlen + 2 + SIMD_WIDTH;
|
||||
int val_mem_size = (col_size * 9 * 2 + 31) >> 5 << 5; // 32字节的整数倍
|
||||
int mem_size = (seq_size + ref_size) * 2 + val_mem_size;
|
||||
|
||||
SIMD_INIT; // 初始化simd用的数据
|
||||
|
||||
assert(init_score > 0);
|
||||
|
||||
// allocate memory
|
||||
mem = thread_mem_request(tmem, mem_size);
|
||||
qtmem = (int16_t *)&mem[0];
|
||||
seq = &qtmem[0];
|
||||
ref = &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 - 1] = target[tlen - 1 - i];
|
||||
}
|
||||
else
|
||||
{
|
||||
for (i = 0; i < qlen; ++i)
|
||||
seq[i] = query[i];
|
||||
for (i = 0; i < tlen; ++i)
|
||||
ref[i + SIMD_WIDTH - 1] = target[i];
|
||||
}
|
||||
|
||||
vmem = &ref[ref_size];
|
||||
for (i = 0; i < (val_mem_size >> 1); i += SIMD_WIDTH)
|
||||
{
|
||||
_mm256_storeu_si256((__m256i *)&vmem[i], zero_vec);
|
||||
}
|
||||
hA = &vmem[0];
|
||||
mA = &vmem[col_size * 3];
|
||||
eA = &vmem[col_size * 5];
|
||||
fA = &vmem[col_size * 7];
|
||||
|
||||
last_max_arr = &hA[0];
|
||||
cur_max_arr = &hA[col_size];
|
||||
next_max_arr = &cur_max_arr[col_size];
|
||||
cur_match_arr = &mA[0];
|
||||
next_match_arr = &mA[col_size];
|
||||
cur_del_arr = &eA[0];
|
||||
next_del_arr = &eA[col_size];
|
||||
cur_ins_arr = &fA[0];
|
||||
next_ins_arr = &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
|
||||
last_max_arr[0] = init_score; // 左上角
|
||||
cur_ins_arr[1] = MAX(0, init_score - (o_ins + e_ins));
|
||||
next_del_arr[0] = init_score;
|
||||
cur_max_arr[1] = cur_ins_arr[1];
|
||||
|
||||
if (qlen == 0 || tlen == 0)
|
||||
Dloop = 0; // 防止意外情况
|
||||
if (window_size >= qlen)
|
||||
{
|
||||
max_ie = 0;
|
||||
gscore = 0;
|
||||
}
|
||||
|
||||
for (D = 1; LIKELY(D < Dloop); ++D)
|
||||
{
|
||||
if (D < tlen)
|
||||
beg1 = 1;
|
||||
else
|
||||
beg1 = D - tlen + 1;
|
||||
if (D < qlen)
|
||||
end1 = D; // 闭区间
|
||||
else
|
||||
end1 = qlen;
|
||||
beg1 = MAX(D - window_size, beg1);
|
||||
end1 = MIN(D + window_size, end1);
|
||||
|
||||
beg = MAX(beg1, beg);
|
||||
end = MIN(end1, end);
|
||||
if (beg > end)
|
||||
break;
|
||||
|
||||
// beg = beg1;
|
||||
// end = end1;
|
||||
|
||||
iend = D - beg; // ref开始计算的位置，倒序
|
||||
span = end - beg;
|
||||
ibeg = iend - span; // 0开始的ref索引位置
|
||||
|
||||
// 每一轮需要记录的数据
|
||||
int m = 0, mj = -1, mi = -1;
|
||||
max_vec = zero_vec;
|
||||
|
||||
// 处理左边界
|
||||
if (beg == 1)
|
||||
{
|
||||
last_max_arr[0] = next_del_arr[0];
|
||||
cur_match_arr[0] = 0;
|
||||
cur_del_arr[0] = MAX(0, init_score - (o_del + e_del * (iend + 1)));
|
||||
}
|
||||
|
||||
for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH)
|
||||
{
|
||||
// 取数据
|
||||
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;
|
||||
}
|
||||
// 处理上边界
|
||||
if (ibeg == 0)
|
||||
{
|
||||
next_ins_arr[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1)));
|
||||
next_max_arr[end + 1] = next_ins_arr[end + 1];
|
||||
next_match_arr[end + 1] = 0;
|
||||
}
|
||||
|
||||
SIMD_FIND_MAX;
|
||||
|
||||
// 注意最后跳出循环j的值
|
||||
j = end + 1;
|
||||
|
||||
if (j == qlen + 1)
|
||||
{
|
||||
max_ie = gscore > next_max_arr[qlen] ? max_ie : ibeg;
|
||||
gscore = gscore > next_max_arr[qlen] ? gscore : next_max_arr[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 = cur_max_arr[j - 1] | next_max_arr[j];
|
||||
if (has_val)
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
beg = j;
|
||||
for (j = end + 1; LIKELY(j >= beg); --j)
|
||||
{
|
||||
int has_val = cur_max_arr[j - 1] | next_max_arr[j];
|
||||
if (has_val)
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
end = j + 1 <= qlen ? j + 1 : qlen;
|
||||
|
||||
// swap m, h, e, f
|
||||
SWAP_DATA_POINTER;
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
|
@ -1,459 +0,0 @@
|
|||
#include <stdlib.h>
|
||||
#include <stdint.h>
|
||||
#include <assert.h>
|
||||
#include <stdio.h>
|
||||
#include <immintrin.h>
|
||||
#include <emmintrin.h>
|
||||
#include "thread_mem.h"
|
||||
#include "common.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
|
||||
|
||||
#define BASE_BYTES 1
|
||||
#define SCORE_BYTES 1
|
||||
#define BOUNDARY_SCORE_NUM 2
|
||||
#define TMP_SCORE_ARRAY_NUM 9
|
||||
#define MEM_ALIGN_BYTES 32
|
||||
#define ALIGN_SHIFT_BITS 5
|
||||
#define SIMD_BYTES 32
|
||||
|
||||
#define AMBIGUOUS_BASE_CODE 4
|
||||
#define AMBIGUOUS_BASE_SCORE -1
|
||||
|
||||
// 32字节对齐（256位）
|
||||
#define align_mem(x) (((x) + 31) >> 5 << 5)
|
||||
#define align_number(x) align_mem(x)
|
||||
|
||||
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] = {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};
|
||||
|
||||
// #define permute_mask _MM_SHUFFLE(0, 1, 2, 3)
|
||||
#define permute_mask 27
|
||||
// 初始化变量
|
||||
#define SIMD_INIT \
|
||||
int oe_del = o_del + e_del, oe_ins = o_ins + e_ins; \
|
||||
__m256i zero_vec; \
|
||||
__m256i max_vec, last_max_vec = _mm256_set1_epi8(init_score); \
|
||||
__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 *)(&read_seq[j])); \
|
||||
__m256i ts_vec = _mm256_loadu_si256((__m256i *)(&ref_seq[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]);
|
||||
|
||||
// cmp_max = _mm256_xor_si256(last_max_vec, cmp_max);
|
||||
// last_max_vec = _mm256_set1_epi8(m);
|
||||
// 找最大值和位置
|
||||
#define SIMD_FIND_MAX \
|
||||
__m256i cmp_max = _mm256_max_epu8(max_vec, last_max_vec); \
|
||||
cmp_max = _mm256_xor_si256(cmp_max, last_max_vec); \
|
||||
cmp_max = _mm256_cmpeq_epi8(cmp_max, zero_vec); \
|
||||
uint32_t cmp_result = _mm256_movemask_epi8(cmp_max); \
|
||||
if (cmp_result != 4294967295) \
|
||||
{ \
|
||||
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]; \
|
||||
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; \
|
||||
for (; mj + 1 < qlen && mi + 1 < tlen; mj++, mi++) \
|
||||
{ \
|
||||
if (read_seq[mj + 2] == ref_seq[mi + 1 + SIMD_WIDTH]) \
|
||||
{ \
|
||||
m += base_match_score; \
|
||||
} \
|
||||
else \
|
||||
{ \
|
||||
break; \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
last_max_vec = _mm256_set1_epi8(m); \
|
||||
}
|
||||
|
||||
// 每轮迭代后，交换数组
|
||||
#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_aligned(thread_mem_t *tmem,
|
||||
int qlen, // query length 待匹配段碱基的query长度
|
||||
const uint8_t *query, // read碱基序列
|
||||
int tlen, // target length reference的长度
|
||||
const uint8_t *target, // reference序列
|
||||
int extend_left, // 是不是向左扩展
|
||||
int o_del, // deletion 错配开始的惩罚系数
|
||||
int e_del, // deletion extension的惩罚系数
|
||||
int o_ins, // insertion 错配开始的惩罚系数
|
||||
int e_ins, // insertion extension的惩罚系数SIMD_BTYES
|
||||
int base_match_score, // 碱基match时的分数
|
||||
int base_mis_score, // 碱基mismatch时的惩罚分数（正数）
|
||||
int window_size, // 提前剪枝系数，w =100 匹配位置和beg的最大距离
|
||||
int end_bonus, // 如果query比对到了最后一个字符，额外奖励分值
|
||||
int init_score, // 该seed的初始得分（完全匹配query的碱基数）
|
||||
int *_qle, // 匹配得到全局最大得分的碱基在query的位置
|
||||
int *_tle, // 匹配得到全局最大得分的碱基在reference的位置
|
||||
int *_gtle, // query全部匹配上的target的长度
|
||||
int *_gscore, // query的端到端匹配得分
|
||||
int *_max_off) // 取得最大得分时在query和reference上位置差的 最大值
|
||||
{
|
||||
uint8_t *mA1, *mA2, *hA0, *hA1, *eA1, *fA1, *hA2, *eA2, *fA2; // hA0保存上上个col的H，其他的保存上个H E F M
|
||||
uint8_t *read_seq, *ref_seq;
|
||||
int i, ibeg, iend, D, j, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
int span, beg1, end1; // 边界条件计算
|
||||
|
||||
uint8_t *mem_addr;
|
||||
int read_size = align_number(qlen * BASE_BYTES + MEM_ALIGN_BYTES);
|
||||
int ref_size = align_number((tlen + SIMD_WIDTH) * BASE_BYTES);
|
||||
int back_diagnal_num = tlen + qlen; // 循环跳出条件 D从1开始遍历
|
||||
int score_array_size = align_number((qlen + BOUNDARY_SCORE_NUM) * SCORE_BYTES);
|
||||
int score_element_num = score_array_size / SCORE_BYTES;
|
||||
int score_mem_size = score_array_size * TMP_SCORE_ARRAY_NUM;
|
||||
int request_mem_size = read_size + ref_size + score_mem_size + MEM_ALIGN_BYTES * 3;
|
||||
|
||||
SIMD_INIT; // 初始化simd用的数据
|
||||
|
||||
assert(init_score > 0);
|
||||
|
||||
mem_addr = thread_mem_request(tmem, request_mem_size);
|
||||
mem_addr = (void *)align_mem((uint64_t)mem_addr);
|
||||
ref_seq = (uint8_t *)&mem_addr[0];
|
||||
read_seq = (uint8_t *)(mem_addr + ref_size + SIMD_BYTES - BASE_BYTES);
|
||||
|
||||
if (extend_left)
|
||||
{
|
||||
for (i = 0; i < qlen; ++i)
|
||||
read_seq[i + 1] = query[qlen - 1 - i];
|
||||
for (i = 0; i < tlen; ++i)
|
||||
ref_seq[i + SIMD_WIDTH - 1] = target[tlen - 1 - i];
|
||||
}
|
||||
else
|
||||
{
|
||||
for (i = 0; i < qlen; ++i)
|
||||
read_seq[i + 1] = query[i];
|
||||
for (i = 0; i < tlen; ++i)
|
||||
ref_seq[i + SIMD_WIDTH - 1] = target[i];
|
||||
}
|
||||
|
||||
mem_addr += read_size + ref_size;
|
||||
for (i = 0; i <= score_mem_size; i += SIMD_BYTES)
|
||||
{
|
||||
_mm256_storeu_si256((__m256i *)&mem_addr[i], zero_vec);
|
||||
}
|
||||
mem_addr += SIMD_BYTES - SCORE_BYTES;
|
||||
|
||||
hA0 = (uint8_t *)&mem_addr[0];
|
||||
hA1 = &hA0[score_element_num];
|
||||
hA2 = &hA1[score_element_num];
|
||||
mA1 = &hA2[score_element_num];
|
||||
mA2 = &mA1[score_element_num];
|
||||
eA1 = &mA2[score_element_num];
|
||||
eA2 = &eA1[score_element_num];
|
||||
fA1 = &eA2[score_element_num];
|
||||
fA2 = &fA1[score_element_num];
|
||||
|
||||
// adjust $window_size if it is too large
|
||||
|
||||
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; // 左上角
|
||||
fA1[1] = MAX(0, init_score - (o_ins + e_ins));
|
||||
eA2[0] = init_score;
|
||||
hA1[1] = fA1[1];
|
||||
|
||||
if (qlen == 0 || tlen == 0)
|
||||
back_diagnal_num = 0; // 防止意外情况
|
||||
if (window_size >= qlen)
|
||||
{
|
||||
max_ie = 0;
|
||||
gscore = 0;
|
||||
}
|
||||
|
||||
for (D = 1; LIKELY(D < back_diagnal_num); ++D)
|
||||
{
|
||||
// 边界条件一定要注意！ tlen 大于，等于，小于 qlen时的情况
|
||||
if (D < tlen)
|
||||
beg1 = 1;
|
||||
else
|
||||
beg1 = D - tlen + 1;
|
||||
if (D < qlen)
|
||||
end1 = D; // 闭区间
|
||||
else
|
||||
end1 = qlen;
|
||||
beg1 = MAX(D - window_size, beg1);
|
||||
end1 = MIN(D + window_size, end1);
|
||||
|
||||
beg = MAX(beg1, beg);
|
||||
end = MIN(end1, end);
|
||||
if (beg > end)
|
||||
break;
|
||||
|
||||
// beg = beg1;
|
||||
// end = end1;
|
||||
|
||||
iend = D - beg; // ref开始计算的位置，倒序
|
||||
span = end - beg;
|
||||
ibeg = iend - span; // 0开始的ref索引位置
|
||||
|
||||
// 每一轮需要记录的数据
|
||||
int m = 0, mj = -1, mi = -1;
|
||||
max_vec = zero_vec;
|
||||
|
||||
// 处理左边界
|
||||
if (beg == 1)
|
||||
{
|
||||
hA0[0] = eA2[0];
|
||||
mA1[0] = 0;
|
||||
eA1[0] = MAX(0, init_score - (o_del + e_del * (iend + 1)));
|
||||
}
|
||||
|
||||
for (j = beg, i = iend; j <= end + 1 - SIMD_WIDTH; j += SIMD_WIDTH, i -= SIMD_WIDTH)
|
||||
{
|
||||
// 取数据
|
||||
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;
|
||||
}
|
||||
|
||||
// 处理上边界
|
||||
if (ibeg == 0)
|
||||
{
|
||||
fA2[end + 1] = MAX(0, init_score - (o_ins + e_ins * (end + 1)));
|
||||
hA2[end + 1] = fA2[end + 1];
|
||||
mA2[end + 1] = 0;
|
||||
}
|
||||
SIMD_FIND_MAX;
|
||||
|
||||
// 注意最后跳出循环j的值
|
||||
j = end + 1;
|
||||
|
||||
if (j == qlen + 1)
|
||||
{
|
||||
max_ie = gscore > hA2[qlen] ? max_ie : ibeg;
|
||||
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;
|
||||
}
|
||||
end = j + 1 <= qlen ? j + 1 : qlen;
|
||||
|
||||
// swap m, h, e, f
|
||||
SWAP_DATA_POINTER;
|
||||
}
|
||||
|
||||
thread_mem_release(tmem, request_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;
|
||||
}
|
||||
|
|
@ -1 +0,0 @@
|
|||
#include "common.h"
|
||||
190
ksw_ext_normal.c
190
ksw_ext_normal.c
|
|
@ -1,190 +0,0 @@
|
|||
#include <stdint.h>
|
||||
#include <stdlib.h>
|
||||
#include <assert.h>
|
||||
#include <stdio.h>
|
||||
#include "thread_mem.h"
|
||||
#include "common.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
|
||||
|
||||
typedef struct
|
||||
{
|
||||
int32_t h, e;
|
||||
} eh_t;
|
||||
|
||||
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;
|
||||
eh_t *eh; // score array
|
||||
int8_t *qp; // query profile
|
||||
int i, j, k, oe_del = o_del + e_del, oe_ins = o_ins + e_ins, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
assert(h0 > 0);
|
||||
// qp = malloc(qlen * m);
|
||||
// eh = calloc(qlen + 1, 8);
|
||||
qp = thread_mem_request(tmem, qlen * m);
|
||||
eh = thread_mem_request_and_clean(tmem, (qlen + 1) * 8);
|
||||
// generate the query profile
|
||||
for (k = i = 0; k < m; ++k)
|
||||
{
|
||||
const int8_t *p = &mat[k * m];
|
||||
for (j = 0; j < qlen; ++j)
|
||||
qp[i++] = p[query[j]]; // 对于qp数组，第0到qlen个元素，表示和A比对的分值，第qlen到2*qlen，表示和C比对的分值，以此类推
|
||||
}
|
||||
// fill the first row
|
||||
// 初始化第一行分值
|
||||
eh[0].h = h0;
|
||||
eh[1].h = h0 > oe_ins ? h0 - oe_ins : 0;
|
||||
for (j = 2; j <= qlen && eh[j - 1].h > e_ins; ++j)
|
||||
eh[j].h = eh[j - 1].h - e_ins;
|
||||
// adjust $w if it is too large
|
||||
k = m * m; // 字符矩阵
|
||||
for (i = 0, max = 0; i < k; ++i) // get the max score
|
||||
max = max > mat[i] ? max : mat[i];
|
||||
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.); // 最大可插入的长度？
|
||||
max_ins = max_ins > 1 ? max_ins : 1;
|
||||
w = w < max_ins ? w : max_ins;
|
||||
max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.);
|
||||
max_del = max_del > 1 ? max_del : 1;
|
||||
w = w < max_del ? w : max_del; // TODO: is this necessary? 上述几行代码都是为了看看能否缩小窗口，减小计算
|
||||
// DP loop
|
||||
max = h0, max_i = max_j = -1;
|
||||
max_ie = -1, gscore = -1;
|
||||
max_off = 0;
|
||||
beg = 0, end = qlen;
|
||||
#ifdef DEBUG_OUT
|
||||
fprintf(score_f_arr[0], "%-4d", h0);
|
||||
fprintf(ins_ext_f_arr[0], "%-4d", h0);
|
||||
fprintf(del_ext_f_arr[0], "%-4d", h0);
|
||||
for (j = 0; LIKELY(j < end); ++j) // 遍历query字符序列
|
||||
{
|
||||
fprintf(score_f_arr[0], "%-4d", MAX(h0 - o_ins - (j + 1) * e_ins, 0));
|
||||
fprintf(ins_ext_f_arr[0], "%-4d", MAX(h0 - o_ins - (j + 1) * e_ins, 0));
|
||||
fprintf(del_ext_f_arr[0], "%-4d", 0);
|
||||
}
|
||||
fprintf(score_f_arr[0], "\n");
|
||||
fprintf(ins_ext_f_arr[0], "\n");
|
||||
fprintf(del_ext_f_arr[0], "\n");
|
||||
#endif
|
||||
for (i = 0; LIKELY(i < tlen); ++i) // 对target逐个字符进行遍历
|
||||
{
|
||||
int t, f = 0, h1, m = 0, mj = -1;
|
||||
// 对于target第i个字符，query中每个字符的分值，只有匹配和不匹配
|
||||
int8_t *q = &qp[target[i] * qlen];
|
||||
// apply the band and the constraint (if provided)
|
||||
if (beg < i - w) // 检查开始点是否可以缩小一些
|
||||
beg = i - w;
|
||||
if (end > i + w + 1) // 检查终点是否可以缩小，使得整体的遍历范围缩小
|
||||
end = i + w + 1;
|
||||
if (end > qlen) // 终点不超过query长度
|
||||
end = qlen;
|
||||
beg = 0;
|
||||
end = qlen;
|
||||
// compute the first column
|
||||
if (beg == 0)
|
||||
{
|
||||
h1 = h0 - (o_del + e_del * (i + 1)); // 只消耗了target序列，query从第一个字符开始匹配第i个target字符
|
||||
if (h1 < 0)
|
||||
h1 = 0;
|
||||
}
|
||||
else
|
||||
h1 = 0;
|
||||
#ifdef DEBUG_OUT
|
||||
fprintf(ins_ext_f_arr[0], "%-4d", 0);
|
||||
fprintf(del_ext_f_arr[0], "%-4d", MAX(h0 - o_del - (i + 1) * e_del, 0));
|
||||
#endif
|
||||
for (j = beg; LIKELY(j < end); ++j) // 遍历query字符序列
|
||||
{
|
||||
#ifdef DEBUG_OUT
|
||||
fprintf(score_f_arr[0], "%-4d", h1);
|
||||
fprintf(ins_ext_f_arr[0], "%-4d", f);
|
||||
fprintf(del_ext_f_arr[0], "%-4d", eh[j].e);
|
||||
#endif
|
||||
// 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 // E表示delete，只消耗target
|
||||
// F(i,j+1) = max{H(i,j)-gapo, F(i,j)} - gape // F表示insert，只消耗query，target的row id固定，query的col index一直增加
|
||||
eh_t *p = &eh[j];
|
||||
int h, M = p->h, e = p->e; // get H(i-1,j-1) and E(i-1,j) // 获取上一轮h值和e值
|
||||
p->h = h1; // set H(i,j-1) for the next row // h1是上一轮计算的结果
|
||||
M = M ? M + q[j] : 0; // separating H and M to disallow a cigar like "100M3I3D20M" // M大于0则当前两个字符进行匹配，无论是否相等，将分值加到M上，此时M可能变为负数
|
||||
h = M > e ? M : e; // e and f are guaranteed to be non-negative, so h>=0 even if M<0 // e和f保证是非负的，所以h肯定非负，即使M可能是负数，因为h取e,f和M的最大值
|
||||
h = h > f ? h : f;
|
||||
h1 = h; // save H(i,j) to h1 for the next column // 用h1来保存当前表格（i,j)对应的分值，用来下次计算
|
||||
mj = m > h ? mj : j; // record the position where max score is achieved // 记录取得最大值时query的字符位置
|
||||
m = m > h ? m : h; // m is stored at eh[mj+1] 因为eh[mj+1]->h表示的是H(i, mj)，及上一轮记录的h
|
||||
t = M - oe_del; // 用来计算delete，假设当前字符(i,j)匹配（无论match还是mismatch），target下一个字符串被空消耗（delete)的分值F(i+1, j)
|
||||
t = t > 0 ? t : 0;
|
||||
e -= e_del; // 假设当前query字符
|
||||
e = e > t ? e : t; // computed E(i+1,j) // t表示(i,j)强行匹配，(i+1, j)是delete的分数，此前e表示(i+1,j)继续delete的分数
|
||||
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)
|
||||
}
|
||||
#ifdef DEBUG_OUT
|
||||
fprintf(score_f_arr[0], "%-4d", h1);
|
||||
fprintf(score_f_arr[0], "\n");
|
||||
fprintf(ins_ext_f_arr[0], "\n");
|
||||
fprintf(del_ext_f_arr[0], "\n");
|
||||
#endif
|
||||
eh[end].h = h1; // end是query序列之外的位置
|
||||
eh[end].e = 0;
|
||||
if (j == qlen) // 此轮遍历到了query的最后一个字符
|
||||
{
|
||||
max_ie = gscore > h1 ? max_ie : i; // max_ie表示取得全局最大分值时，target字符串的位置
|
||||
gscore = gscore > h1 ? gscore : h1;
|
||||
}
|
||||
// if (m == 0) // 遍历完query之后，当前轮次的最大分值为0，则跳出循环
|
||||
// break;
|
||||
if (m > max) // 当前轮最大分值大于之前的最大分值
|
||||
{
|
||||
max = m, max_i = i, max_j = mj; // 更新取得最大值的target和query的位置
|
||||
max_off = max_off > abs(mj - i) ? max_off : abs(mj - i); // 取得最大分值时，query和target对应字符串坐标的差值
|
||||
}
|
||||
else if (0) //(zdrop > 0) // 当前轮匹配之后取得的最大分值没有大于之前的最大值，而且zdrop值大于0
|
||||
{
|
||||
if (i - max_i > mj - max_j)
|
||||
{
|
||||
if (max - m - ((i - max_i) - (mj - max_j)) * e_del > zdrop) // 之前最大分值 -（从取得最大值的点出发，当前的delete总长度对应的分值 + 当前轮取得的最大值） > zdrop
|
||||
break;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (max - m - ((mj - max_j) - (i - max_i)) * e_ins > zdrop) // 同上，不过这次是insert，可能是说明有很多mismatch
|
||||
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; // 剪枝没考虑f，即insert
|
||||
// beg = 0, end = qlen; // uncomment this line for debugging
|
||||
// fprintf(stderr, "\n");
|
||||
// fprintf(stderr, "%d\n", end);
|
||||
}
|
||||
// free(eh);
|
||||
// free(qp);
|
||||
thread_mem_release(tmem, qlen * m + (qlen + 1) * 8);
|
||||
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;
|
||||
}
|
||||
129
main.c
129
main.c
|
|
@ -14,7 +14,7 @@
|
|||
#include <time.h>
|
||||
#include "sys/time.h"
|
||||
#include "thread_mem.h"
|
||||
#include "ksw_ext.h"
|
||||
#include "bsw.h"
|
||||
#include "utils.h"
|
||||
#include "common.h"
|
||||
|
||||
|
|
@ -38,47 +38,8 @@ int64_t get_mseconds()
|
|||
int64_t time_sw[KERNEL_NUM] = {0};
|
||||
#endif
|
||||
|
||||
#ifdef DEBUG_RETURN_VALUE
|
||||
#define OUTPUT_RETVAL_1(kernel_num) \
|
||||
fprintf(retval_f_arr[kernel_num], "%d\t%d\t%d\t%d\t", score[kernel_num], info_arr[i][0], info_arr[i][1], info_arr[i][2]); \
|
||||
for (j = cur_query_pos - info_arr[i][0]; j < cur_query_pos; ++j) \
|
||||
{ \
|
||||
fprintf(retval_f_arr[kernel_num], "%c", t_2bit2char[(uint8_t)query_arr[j]]); \
|
||||
} \
|
||||
fprintf(retval_f_arr[kernel_num], "\t"); \
|
||||
for (j = cur_target_pos - info_arr[i][1]; j < cur_target_pos; ++j) \
|
||||
{ \
|
||||
fprintf(retval_f_arr[kernel_num], "%c", t_2bit2char[(uint8_t)target_arr[j]]); \
|
||||
} \
|
||||
fprintf(retval_f_arr[kernel_num], "\n");
|
||||
|
||||
#define OUTPUT_RETVAL_SCORE \
|
||||
fprintf(retval_f_arr[kernel_num], "%d\t%d\t%d\t%d\t%d\t%d\n", score[kernel_num], qle, tle, gtle, gscore, max_off[0])
|
||||
// fprintf(retval_f_arr[kernel_num], " %d %d\n", cur_query_pos, info_arr[i][0]);
|
||||
// fprintf(retval_f_arr[kernel_num], "%d\t%d\t%d\t%d\n", score[kernel_num], info_arr[i][0], info_arr[i][1], info_arr[i][2])
|
||||
|
||||
#define OUTPUT_RETVAL_target(kernel_num) \
|
||||
for (j = cur_target_pos - info_arr[i][1]; j < cur_target_pos; ++j) \
|
||||
{ \
|
||||
fprintf(retval_f_arr[kernel_num], "%c", t_2bit2char[(uint8_t)target_arr[j]]); \
|
||||
} \
|
||||
fprintf(retval_f_arr[kernel_num], "\n");
|
||||
|
||||
#define OUTPUT_RETVAL_query(kernel_num) \
|
||||
for (j = cur_query_pos - info_arr[i][0]; j < cur_query_pos; ++j) \
|
||||
{ \
|
||||
fprintf(retval_f_arr[kernel_num], "%c", t_2bit2char[(uint8_t)query_arr[j]]); \
|
||||
} \
|
||||
fprintf(retval_f_arr[kernel_num], "\n");
|
||||
|
||||
#define OUTPUT_RETVAL_INFO(kernel_num) \
|
||||
fprintf(retval_f_arr[kernel_num], "%-8d%-8d%-8d\n", info_arr[i][0], info_arr[i][1], info_arr[i][2]);
|
||||
#define OUTPUT_RETVAL(kernel_num) OUTPUT_RETVAL_target(kernel_num)
|
||||
#else
|
||||
#define OUTPUT_RETVAL(kernel_num)
|
||||
#endif
|
||||
|
||||
#define _PERFORMANCE_TEST_NORMAL(kernel_num, func) \
|
||||
#define _PERFORMANCE_TEST(kernel_num, func) \
|
||||
cur_query_pos = 0; \
|
||||
cur_target_pos = 0; \
|
||||
for (i = 0; i < block_line_num; ++i) \
|
||||
|
|
@ -95,52 +56,16 @@ int64_t time_sw[KERNEL_NUM] = {0};
|
|||
score_total[kernel_num] += score[kernel_num]; \
|
||||
cur_query_pos += info_arr[i][0]; \
|
||||
cur_target_pos += info_arr[i][1]; \
|
||||
OUTPUT_RETVAL(0); \
|
||||
}
|
||||
|
||||
#define _PERFORMANCE_TEST_AVX2(kernel_num, func) \
|
||||
cur_query_pos = 0; \
|
||||
cur_target_pos = 0; \
|
||||
for (i = 0; i < block_line_num; ++i) \
|
||||
{ \
|
||||
score[kernel_num] = func( \
|
||||
&tmem[kernel_num], \
|
||||
info_arr[i][0], \
|
||||
(uint8_t *)query_arr + cur_query_pos, \
|
||||
info_arr[i][1], \
|
||||
(uint8_t *)target_arr + cur_target_pos, \
|
||||
0, 6, 1, 6, 1, \
|
||||
1, 4, \
|
||||
100, 5, \
|
||||
info_arr[i][2], \
|
||||
&qle, &tle, >le, &gscore, &max_off[0]); \
|
||||
score_total[kernel_num] += score[kernel_num]; \
|
||||
cur_query_pos += info_arr[i][0]; \
|
||||
cur_target_pos += info_arr[i][1]; \
|
||||
OUTPUT_RETVAL(kernel_num); \
|
||||
}
|
||||
#define _PERFORMANCE_TEST_AVX2_T1(kernel_num, func) \
|
||||
cur_query_pos = 0; \
|
||||
cur_target_pos = 0; \
|
||||
for (i = 0; i < block_line_num; ++i) \
|
||||
{ \
|
||||
cur_query_pos += info_arr[i][0]; \
|
||||
cur_target_pos += info_arr[i][1]; \
|
||||
OUTPUT_RETVAL(kernel_num); \
|
||||
}
|
||||
|
||||
#ifdef SHOW_PERF
|
||||
#define PERFORMANCE_TEST_NORMAL(kernel_num, func) \
|
||||
start_time = get_mseconds(); \
|
||||
_PERFORMANCE_TEST_NORMAL(kernel_num, func); \
|
||||
time_sw[kernel_num] += get_mseconds() - start_time
|
||||
|
||||
#define PERFORMANCE_TEST_AVX2(kernel_num, func) \
|
||||
start_time = get_mseconds(); \
|
||||
_PERFORMANCE_TEST_AVX2_T1(kernel_num, func); \
|
||||
#define PERFORMANCE_TEST(kernel_num, func) \
|
||||
start_time = get_mseconds(); \
|
||||
_PERFORMANCE_TEST(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)
|
||||
#define PERFORMANCE_TEST(kernel_num, func) _PERFORMANCE_TEST(kernel_num, func)
|
||||
#endif
|
||||
|
||||
// 读取一行序列数据
|
||||
|
|
@ -245,6 +170,7 @@ int main(int argc, char *argv[])
|
|||
j += 3;
|
||||
}
|
||||
|
||||
int64_t all_qlen = 0;
|
||||
while (!feof(target_f))
|
||||
{
|
||||
block_line_num = 0; // 记录每次读取的行数
|
||||
|
|
@ -284,6 +210,7 @@ int main(int argc, char *argv[])
|
|||
}
|
||||
|
||||
// 读info
|
||||
|
||||
cur_read_size = 0;
|
||||
for (i = 0; i < block_line_num; ++i)
|
||||
{
|
||||
|
|
@ -295,54 +222,52 @@ int main(int argc, char *argv[])
|
|||
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);
|
||||
all_qlen += info_arr[i][0];
|
||||
}
|
||||
|
||||
#ifdef DEBUG_RETURN_VALUE
|
||||
cur_read_size = 0;
|
||||
for (i = 0; i < block_line_num; ++i)
|
||||
{
|
||||
fprintf(stdout, "%d\n", i);
|
||||
for (j = cur_read_size; j < cur_read_size + info_arr[i][1]; ++j)
|
||||
{
|
||||
fprintf(stdout, "j: %d\n", j);
|
||||
fprintf(retval_f_arr[0], "%c", t_2bit2char[(uint8_t)target_arr[j]]);
|
||||
}
|
||||
fprintf(retval_f_arr[0], "\n");
|
||||
cur_read_size += info_arr[i][1];
|
||||
}
|
||||
#endif
|
||||
// for (i = 0; i < block_line_num; ++i)
|
||||
//{
|
||||
// fprintf(retval_f_arr[0], "%d\n", info_arr[i][1]);
|
||||
//}
|
||||
|
||||
// 性能测试
|
||||
|
||||
#if 1
|
||||
// normal sw
|
||||
// PERFORMANCE_TEST_NORMAL(0, ksw_extend_normal);
|
||||
PERFORMANCE_TEST(0, normal);
|
||||
|
||||
// normal pruning
|
||||
PERFORMANCE_TEST(1, normal_pruning);
|
||||
|
||||
// avx2
|
||||
// PERFORMANCE_TEST_AVX2(1, ksw_extend_avx2);
|
||||
// avx2 heuristics
|
||||
// PERFORMANCE_TEST_AVX2(2, ksw_extend_avx2_heuristics);
|
||||
// avx2 mem aligned
|
||||
// PERFORMANCE_TEST_AVX2(3, ksw_extend_avx2_aligned);
|
||||
|
||||
// avx2 u8
|
||||
// PERFORMANCE_TEST_AVX2(4, ksw_extend_avx2_u8);
|
||||
// avx2 u8 heuristics
|
||||
// PERFORMANCE_TEST_AVX2(5, ksw_extend_avx2_u8_heuristics);
|
||||
// avx2 u8 mem aligned
|
||||
// PERFORMANCE_TEST_AVX2(6, ksw_extend_avx2_u8_aligned);
|
||||
PERFORMANCE_TEST(2, avx2_u8);
|
||||
// avx2 pruning
|
||||
PERFORMANCE_TEST(3, avx2_u8_pruning);
|
||||
#endif
|
||||
}
|
||||
|
||||
fprintf(stderr, "%ld\n", total_line_num);
|
||||
fprintf(stderr, "all_qlen: %ld\n", all_qlen);
|
||||
|
||||
#ifdef SHOW_PERF
|
||||
char *kernel_names[7] = {
|
||||
char *kernel_names[4] = {
|
||||
"normal",
|
||||
"avx2",
|
||||
"avx2_heuristics",
|
||||
"avx2_aligned",
|
||||
"normal_pruning",
|
||||
"avx2_u8",
|
||||
"avx2_u8_heuristics",
|
||||
"avx2_u8_aligned"};
|
||||
"avx2_u8_pruning"};
|
||||
|
||||
for (i = 0; i < KERNEL_NUM; ++i)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -0,0 +1,162 @@
|
|||
#include <stdint.h>
|
||||
#include <stdlib.h>
|
||||
#include <assert.h>
|
||||
#include <stdio.h>
|
||||
#include "thread_mem.h"
|
||||
#include "common.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
|
||||
|
||||
typedef struct
|
||||
{
|
||||
int32_t h, e;
|
||||
} eh_t;
|
||||
|
||||
int 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;
|
||||
eh_t *eh; // score array
|
||||
int8_t *qp; // query profile
|
||||
int i, j, k, oe_del = o_del + e_del, oe_ins = o_ins + e_ins, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
assert(h0 > 0);
|
||||
qp = malloc(qlen * m);
|
||||
eh = calloc(qlen + 1, 8);
|
||||
//qp = thread_mem_request(tmem, qlen * m);
|
||||
//eh = thread_mem_request_and_clean(tmem, (qlen + 1) * 8);
|
||||
// generate the query profile
|
||||
//fprintf(stdout, "%d\n", h0);
|
||||
|
||||
for (k = i = 0; k < m; ++k)
|
||||
{
|
||||
const int8_t *p = &mat[k * m];
|
||||
for (j = 0; j < qlen; ++j)
|
||||
qp[i++] = p[query[j]]; // 对于qp数组，第0到qlen个元素，表示和A比对的分值，第qlen到2*qlen，表示和C比对的分值，以此类推
|
||||
}
|
||||
// fill the first row
|
||||
// 初始化第一行分值
|
||||
eh[0].h = h0;
|
||||
eh[1].h = h0 > oe_ins ? h0 - oe_ins : 0;
|
||||
for (j = 2; j <= qlen && eh[j - 1].h > e_ins; ++j)
|
||||
eh[j].h = eh[j - 1].h - e_ins;
|
||||
// adjust $w if it is too large
|
||||
k = m * m; // 字符矩阵
|
||||
for (i = 0, max = 0; i < k; ++i) // get the max score
|
||||
max = max > mat[i] ? max : mat[i];
|
||||
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.); // 最大可插入的长度？
|
||||
max_ins = max_ins > 1 ? max_ins : 1;
|
||||
w = w < max_ins ? w : max_ins;
|
||||
max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.);
|
||||
max_del = max_del > 1 ? max_del : 1;
|
||||
w = w < max_del ? w : max_del; // TODO: is this necessary? 上述几行代码都是为了看看能否缩小窗口，减小计算
|
||||
// DP loop
|
||||
max = h0, max_i = max_j = -1;
|
||||
max_ie = -1, gscore = -1;
|
||||
max_off = 0;
|
||||
beg = 0, end = qlen;
|
||||
for (i = 0; LIKELY(i < tlen); ++i) // 对target逐个字符进行遍历
|
||||
{
|
||||
int t, f = 0, h1, m = 0, mj = -1;
|
||||
// 对于target第i个字符，query中每个字符的分值，只有匹配和不匹配
|
||||
int8_t *q = &qp[target[i] * qlen];
|
||||
// apply the band and the constraint (if provided)
|
||||
if (beg < i - w) // 检查开始点是否可以缩小一些
|
||||
beg = i - w;
|
||||
if (end > i + w + 1) // 检查终点是否可以缩小，使得整体的遍历范围缩小
|
||||
end = i + w + 1;
|
||||
if (end > qlen) // 终点不超过query长度
|
||||
end = qlen;
|
||||
// compute the first column
|
||||
if (beg == 0)
|
||||
{
|
||||
h1 = h0 - (o_del + e_del * (i + 1)); // 只消耗了target序列，query从第一个字符开始匹配第i个target字符
|
||||
if (h1 < 0)
|
||||
h1 = 0;
|
||||
}
|
||||
else
|
||||
h1 = 0;
|
||||
|
||||
for (j = beg; LIKELY(j < end); ++j) // 遍历query字符序列
|
||||
{
|
||||
// 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 // E表示delete，只消耗target
|
||||
// F(i,j+1) = max{H(i,j)-gapo, F(i,j)} - gape // F表示insert，只消耗query，target的row id固定，query的col index一直增加
|
||||
eh_t *p = &eh[j];
|
||||
int h, M = p->h, e = p->e; // get H(i-1,j-1) and E(i-1,j) // 获取上一轮h值和e值
|
||||
p->h = h1; // set H(i,j-1) for the next row // h1是上一轮计算的结果
|
||||
M = M ? M + q[j] : 0; // separating H and M to disallow a cigar like "100M3I3D20M" // M大于0则当前两个字符进行匹配，无论是否相等，将分值加到M上，此时M可能变为负数
|
||||
h = M > e ? M : e; // e and f are guaranteed to be non-negative, so h>=0 even if M<0 // e和f保证是非负的，所以h肯定非负，即使M可能是负数，因为h取e,f和M的最大值
|
||||
h = h > f ? h : f;
|
||||
h1 = h; // save H(i,j) to h1 for the next column // 用h1来保存当前表格（i,j)对应的分值，用来下次计算
|
||||
mj = m > h ? mj : j; // record the position where max score is achieved // 记录取得最大值时query的字符位置
|
||||
m = m > h ? m : h; // m is stored at eh[mj+1] 因为eh[mj+1]->h表示的是H(i, mj)，及上一轮记录的h
|
||||
t = M - oe_del; // 用来计算delete，假设当前字符(i,j)匹配（无论match还是mismatch），target下一个字符串被空消耗（delete)的分值F(i+1, j)
|
||||
t = t > 0 ? t : 0;
|
||||
e -= e_del; // 假设当前query字符
|
||||
e = e > t ? e : t; // computed E(i+1,j) // t表示(i,j)强行匹配，(i+1, j)是delete的分数，此前e表示(i+1,j)继续delete的分数
|
||||
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; // end是query序列之外的位置
|
||||
eh[end].e = 0;
|
||||
if (j == qlen) // 此轮遍历到了query的最后一个字符
|
||||
{
|
||||
max_ie = gscore > h1 ? max_ie : i; // max_ie表示取得全局最大分值时，target字符串的位置
|
||||
gscore = gscore > h1 ? gscore : h1;
|
||||
}
|
||||
//if (m == 0) // 遍历完query之后，当前轮次的最大分值为0，则跳出循环
|
||||
// break;
|
||||
if (m > max) // 当前轮最大分值大于之前的最大分值
|
||||
{
|
||||
max = m, max_i = i, max_j = mj; // 更新取得最大值的target和query的位置
|
||||
max_off = max_off > abs(mj - i) ? max_off : abs(mj - i); // 取得最大分值时，query和target对应字符串坐标的差值
|
||||
}
|
||||
else if (zdrop > 0) // 当前轮匹配之后取得的最大分值没有大于之前的最大值，而且zdrop值大于0
|
||||
{
|
||||
#if 0
|
||||
if (i - max_i > mj - max_j)
|
||||
{
|
||||
if (max - m - ((i - max_i) - (mj - max_j)) * e_del > zdrop) // 之前最大分值 -（从取得最大值的点出发，当前的delete总长度对应的分值 + 当前轮取得的最大值） > zdrop
|
||||
break;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (max - m - ((mj - max_j) - (i - max_i)) * e_ins > zdrop) // 同上，不过这次是insert，可能是说明有很多mismatch
|
||||
break;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
// update beg and end for the next round
|
||||
//for (j = beg; LIKELY(j < end) && eh[j].h == 0 && eh[j].e == 0; ++j)
|
||||
// ;
|
||||
//beg = j;
|
||||
//for (j = end; LIKELY(j >= beg) && eh[j].h == 0 && eh[j].e == 0; --j)
|
||||
// ;
|
||||
//end = j + 2 < qlen ? j + 2 : qlen; // 剪枝没考虑f，即insert
|
||||
beg = 0, end = qlen; // uncomment this line for debugging
|
||||
}
|
||||
free(eh);
|
||||
free(qp);
|
||||
// thread_mem_release(tmem, qlen * m + (qlen + 1) * 8);
|
||||
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;
|
||||
}
|
||||
|
|
@ -0,0 +1,160 @@
|
|||
#include <stdint.h>
|
||||
#include <stdlib.h>
|
||||
#include <assert.h>
|
||||
#include <stdio.h>
|
||||
#include "thread_mem.h"
|
||||
#include "common.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
|
||||
|
||||
typedef struct
|
||||
{
|
||||
int32_t h, e;
|
||||
} eh_t;
|
||||
|
||||
int normal_pruning(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;
|
||||
eh_t *eh; // score array
|
||||
int8_t *qp; // query profile
|
||||
int i, j, k, oe_del = o_del + e_del, oe_ins = o_ins + e_ins, beg, end, max, max_i, max_j, max_ins, max_del, max_ie, gscore, max_off;
|
||||
assert(h0 > 0);
|
||||
// qp = malloc(qlen * m);
|
||||
// eh = calloc(qlen + 1, 8);
|
||||
qp = thread_mem_request(tmem, qlen * m);
|
||||
eh = thread_mem_request_and_clean(tmem, (qlen + 1) * 8);
|
||||
// generate the query profile
|
||||
for (k = i = 0; k < m; ++k)
|
||||
{
|
||||
const int8_t *p = &mat[k * m];
|
||||
for (j = 0; j < qlen; ++j)
|
||||
qp[i++] = p[query[j]]; // 对于qp数组，第0到qlen个元素，表示和A比对的分值，第qlen到2*qlen，表示和C比对的分值，以此类推
|
||||
}
|
||||
// fill the first row
|
||||
// 初始化第一行分值
|
||||
eh[0].h = h0;
|
||||
eh[1].h = h0 > oe_ins ? h0 - oe_ins : 0;
|
||||
for (j = 2; j <= qlen && eh[j - 1].h > e_ins; ++j)
|
||||
eh[j].h = eh[j - 1].h - e_ins;
|
||||
// adjust $w if it is too large
|
||||
k = m * m; // 字符矩阵
|
||||
for (i = 0, max = 0; i < k; ++i) // get the max score
|
||||
max = max > mat[i] ? max : mat[i];
|
||||
max_ins = (int)((double)(qlen * max + end_bonus - o_ins) / e_ins + 1.); // 最大可插入的长度？
|
||||
max_ins = max_ins > 1 ? max_ins : 1;
|
||||
w = w < max_ins ? w : max_ins;
|
||||
max_del = (int)((double)(qlen * max + end_bonus - o_del) / e_del + 1.);
|
||||
max_del = max_del > 1 ? max_del : 1;
|
||||
w = w < max_del ? w : max_del; // TODO: is this necessary? 上述几行代码都是为了看看能否缩小窗口，减小计算
|
||||
// DP loop
|
||||
max = h0, max_i = max_j = -1;
|
||||
max_ie = -1, gscore = -1;
|
||||
max_off = 0;
|
||||
beg = 0, end = qlen;
|
||||
for (i = 0; LIKELY(i < tlen); ++i) // 对target逐个字符进行遍历
|
||||
{
|
||||
int t, f = 0, h1, m = 0, mj = -1;
|
||||
// 对于target第i个字符，query中每个字符的分值，只有匹配和不匹配
|
||||
int8_t *q = &qp[target[i] * qlen];
|
||||
// apply the band and the constraint (if provided)
|
||||
if (beg < i - w) // 检查开始点是否可以缩小一些
|
||||
beg = i - w;
|
||||
if (end > i + w + 1) // 检查终点是否可以缩小，使得整体的遍历范围缩小
|
||||
end = i + w + 1;
|
||||
if (end > qlen) // 终点不超过query长度
|
||||
end = qlen;
|
||||
// compute the first column
|
||||
if (beg == 0)
|
||||
{
|
||||
h1 = h0 - (o_del + e_del * (i + 1)); // 只消耗了target序列，query从第一个字符开始匹配第i个target字符
|
||||
if (h1 < 0)
|
||||
h1 = 0;
|
||||
}
|
||||
else
|
||||
h1 = 0;
|
||||
|
||||
for (j = beg; LIKELY(j < end); ++j) // 遍历query字符序列
|
||||
{
|
||||
// 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 // E表示delete，只消耗target
|
||||
// F(i,j+1) = max{H(i,j)-gapo, F(i,j)} - gape // F表示insert，只消耗query，target的row id固定，query的col index一直增加
|
||||
eh_t *p = &eh[j];
|
||||
int h, M = p->h, e = p->e; // get H(i-1,j-1) and E(i-1,j) // 获取上一轮h值和e值
|
||||
p->h = h1; // set H(i,j-1) for the next row // h1是上一轮计算的结果
|
||||
M = M ? M + q[j] : 0; // separating H and M to disallow a cigar like "100M3I3D20M" // M大于0则当前两个字符进行匹配，无论是否相等，将分值加到M上，此时M可能变为负数
|
||||
h = M > e ? M : e; // e and f are guaranteed to be non-negative, so h>=0 even if M<0 // e和f保证是非负的，所以h肯定非负，即使M可能是负数，因为h取e,f和M的最大值
|
||||
h = h > f ? h : f;
|
||||
h1 = h; // save H(i,j) to h1 for the next column // 用h1来保存当前表格（i,j)对应的分值，用来下次计算
|
||||
mj = m > h ? mj : j; // record the position where max score is achieved // 记录取得最大值时query的字符位置
|
||||
m = m > h ? m : h; // m is stored at eh[mj+1] 因为eh[mj+1]->h表示的是H(i, mj)，及上一轮记录的h
|
||||
t = M - oe_del; // 用来计算delete，假设当前字符(i,j)匹配（无论match还是mismatch），target下一个字符串被空消耗（delete)的分值F(i+1, j)
|
||||
t = t > 0 ? t : 0;
|
||||
e -= e_del; // 假设当前query字符
|
||||
e = e > t ? e : t; // computed E(i+1,j) // t表示(i,j)强行匹配，(i+1, j)是delete的分数，此前e表示(i+1,j)继续delete的分数
|
||||
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; // end是query序列之外的位置
|
||||
eh[end].e = 0;
|
||||
if (j == qlen) // 此轮遍历到了query的最后一个字符
|
||||
{
|
||||
max_ie = gscore > h1 ? max_ie : i; // max_ie表示取得全局最大分值时，target字符串的位置
|
||||
gscore = gscore > h1 ? gscore : h1;
|
||||
}
|
||||
if (m == 0) // 遍历完query之后，当前轮次的最大分值为0，则跳出循环
|
||||
break;
|
||||
if (m > max) // 当前轮最大分值大于之前的最大分值
|
||||
{
|
||||
max = m, max_i = i, max_j = mj; // 更新取得最大值的target和query的位置
|
||||
max_off = max_off > abs(mj - i) ? max_off : abs(mj - i); // 取得最大分值时，query和target对应字符串坐标的差值
|
||||
}
|
||||
else if (zdrop > 0) // 当前轮匹配之后取得的最大分值没有大于之前的最大值，而且zdrop值大于0
|
||||
{
|
||||
#if 0
|
||||
if (i - max_i > mj - max_j)
|
||||
{
|
||||
if (max - m - ((i - max_i) - (mj - max_j)) * e_del > zdrop) // 之前最大分值 -（从取得最大值的点出发，当前的delete总长度对应的分值 + 当前轮取得的最大值） > zdrop
|
||||
break;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (max - m - ((mj - max_j) - (i - max_i)) * e_ins > zdrop) // 同上，不过这次是insert，可能是说明有很多mismatch
|
||||
break;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
// update beg and end for the next round
|
||||
for (j = beg; LIKELY(j < end) && eh[j].h == 0 && eh[j].e == 0; ++j)
|
||||
;
|
||||
beg = j;
|
||||
for (j = end; LIKELY(j >= beg) && eh[j].h == 0 && eh[j].e == 0; --j)
|
||||
;
|
||||
end = j + 2 < qlen ? j + 2 : qlen; // 剪枝没考虑f，即insert
|
||||
// beg = 0, end = qlen; // uncomment this line for debugging
|
||||
}
|
||||
// free(eh);
|
||||
// free(qp);
|
||||
thread_mem_release(tmem, qlen * m + (qlen + 1) * 8);
|
||||
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;
|
||||
}
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
#!/bin/bash
|
||||
#/home/zzh/work/sw_perf/sw_perf /home/zzh/work/sw_perf/ext_out/q_0.fa /home/zzh/work/sw_perf/ext_out/t_0.fa /home/zzh/work/sw_perf/ext_out/i_0.txt
|
||||
#/home/zzh/work/sw_perf/sw_perf /home/zzh/work/sw_perf/ext_out/q_1.fa /home/zzh/work/sw_perf/ext_out/t_1.fa /home/zzh/work/sw_perf/ext_out/i_1.txt
|
||||
/home/zzh/work/sw_perf/sw_perf /home/zzh/work/sw_perf/ext_out/q_2.fa /home/zzh/work/sw_perf/ext_out/t_2.fa /home/zzh/work/sw_perf/ext_out/i_2.txt
|
||||
#/home/zzh/work/sw_perf/sw_perf /home/zzh/work/sw_perf/ext_out/q_3.fa /home/zzh/work/sw_perf/ext_out/t_3.fa /home/zzh/work/sw_perf/ext_out/i_3.txt
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
#!/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
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
#!/bin/bash
|
||||
#/home/zzh/work/sw_perf/sw_perf /home/zzh/work/sw_perf/input/q.fa /home/zzh/work/sw_perf/input/t.fa /home/zzh/work/sw_perf/input/i.txt
|
||||
/home/zzh/work/sw_perf/sw_perf /home/zzh/work/sw_perf/input/q_d.fa /home/zzh/work/sw_perf/input/t_d.fa /home/zzh/work/sw_perf/input/i_d.txt
|
||||
|
||||
2
run_l.sh
2
run_l.sh
|
|
@ -1,2 +0,0 @@
|
|||
#!/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
|
||||
2
run_m.sh
2
run_m.sh
|
|
@ -1,2 +0,0 @@
|
|||
#!/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
|
||||
Loading…
Reference in New Issue