學校課程/平行程式

批改娘 10112. Longest Common Subsequence (CUDA)

contents

  1. 1. 題目描述
  2. 2. 輸入格式
  3. 3. 輸出格式
  4. 4. 範例輸入
  5. 5. 範例輸出
  6. 6. 編譯參數
  7. 7. Solution

題目描述

給兩個字串 $X, \; Y$，在兩個字串中都有出現且最長的子序列 (subsequence)，就是最長共同子字串
。

輸入格式

有多組測資，每組測資有兩行字串 $X, \; Y$，$X, \; Y$ 只由 A T C G 四個字母構成。

  • $1 \le |X|, |Y| \le 60000$

輸出格式

針對每一組測資，輸出一行 $X, \; Y$ 的最長共同子字串長度。

範例輸入

1
2
3
4
TCA
GTA
TGGAC
TATCT

範例輸出

1
2
2
3

編譯參數

1
$ nvcc -Xcompiler "-O2 -fopenmp" main.cu -o main

Solution

相比於 OpenMP 版本，由於 GPU 的快取設計不如 CPU，在我們使用失敗指針優化 branch 的發生，對 GPU 反而還是慢的，因為要存取 global memory，這使得還不如發生 branch，讓 warp 分多次跑反而比較快，因為有足夠多的 warp 在一個 block。在優化程度上，CPU 和 GPU 還是得分開解決。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include <cstdio>
#include <cstdlib>
#define MAXN 60005
#define CHARSET 4
typedef unsigned short uint16;
static char A[MAXN], B[MAXN];
static int c2i[128];
static char i2c[CHARSET+1] = "ACGT";
static int *cuP;
static char *cuB;
static uint16 *cuDP;
__global__ void prepare(int *P, char *B, int nb) {
    int i = threadIdx.x;
    int *p = P + i*MAXN;
    char c = "ACGT"[i];
    p[0] = 0;
    for (int j = 1; j <= nb; j++)
        p[j] = (B[j] == c) ? j : p[j-1];
}
__global__ void run(int nb, uint16 *dpIn, uint16 *dpOut, int *P) {
    int j = blockDim.x*blockIdx.x+threadIdx.x+1;
    if (j > nb)	return ;
    int pos = P[j];
    uint16 t1 = pos ? dpIn[pos-1]+1 : 0;
    uint16 t2 = dpIn[j];
    dpOut[j] = t1 > t2 ? t1 : t2;
}
int lcs_len(const char *A, int na, const char *B, int nb, uint16 dpf[]) {
    B--;
    cudaMemcpy(cuB, B, sizeof(char)*(nb+1), cudaMemcpyHostToDevice);
    cudaMemset(cuDP, 0, sizeof(uint16)*(nb+1));
    cudaMemset(cuDP+MAXN, 0, sizeof(uint16)*(nb+1));
    
    int bsz = 1024;
    dim3 bb(bsz);
    dim3 gg((nb+bsz-1)/bsz);
    prepare<<<1, 4>>>(cuP, cuB, nb);
    for (int i = 0; i < na; i++) {
        int v = c2i[A[i]];
        run<<<gg, bb>>>(nb, cuDP+(i&1)*MAXN, cuDP+((i&1)^1)*MAXN, cuP+v*MAXN);
    }
    cudaMemcpy(dpf, cuDP+(na&1)*MAXN, sizeof(uint16)*(nb+1), cudaMemcpyDeviceToHost);
    return dpf[nb];
}
int main() {
    for (int i = 0; i < CHARSET; i++)
        c2i[i2c[i]] = i;
    cudaMalloc(&cuB, sizeof(char)*MAXN);
    cudaMalloc(&cuP, sizeof(int)*MAXN*4);
    cudaMalloc(&cuDP, sizeof(uint16)*MAXN*2);
    static uint16 dpf[MAXN];
    while (scanf("%s %s", A, B) == 2) {
        int len = lcs_len(A, strlen(A), B, strlen(B), dpf);
        printf("%d\n", len);
    }
    cudaFree(cuB);
    cudaFree(cuP);
    cudaFree(cuDP);
    return 0;
}