2016-06-28

批改娘 10101. Fast Game of Life (CUDA)

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

題目描述

生命遊戲中，對於任意細胞，規則如下：
每個細胞有兩種狀態-存活或死亡，每個細胞與以自身為中心的周圍八格細胞產生互動。

當前細胞為存活狀態時，當周圍低於 2 個 (不包含 2 個) 存活細胞時，該細胞變成死亡狀態。
當前細胞為存活狀態時，當周圍有 2 個或 3 個存活細胞時，該細胞保持原樣。
當前細胞為存活狀態時，當周圍有 3 個以上的存活細胞時，該細胞變成死亡狀態。
當前細胞為死亡狀態時，當周圍有 3 個存活細胞時，該細胞變成存活狀態。

可以把最初的細胞結構定義為種子，當所有在種子中的細胞同時被以上規則處理後，可以得到第一代細胞圖。按規則繼續處理當前的細胞圖，可以得到下一代的細胞圖，周而復始。

輸入格式

輸入第一行有兩個整數 $N$, $M$，表示盤面大小為 $N \times N$，模擬週期次數 $M$。接下來會有 $N$ 行，每一行上會有 $N$ 個字符，以 0 表示 $(i, j)$ 格子上的細胞屬於死亡狀態，反之 1 為存活狀態。

$1 \le N \le 2000$
$0 \le M \le 5000$

輸出格式

對於每一組測資輸出 $N$ 行，每一行上有 $N$ 個字元表示模擬 $M$ 次的最終盤面結果。

範例輸入 1

範例輸出 1

範例輸入 2

範例輸出 2

編譯參數

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

Solution

CUDA 的函數分成 blocking 或者是 non-blocking，這是因為是異質計算在 host 上不一定要等到 GPU 算完才能執行下一行指令。kernel function call 是 non-blocking 的，可以藉由 cudaDeviceSynchronize 等待 device 所有的 task 都完成，才進行到下一個運行區塊。

但是別忘記 cudaMemcpy 和 kernel function call 這一類都類似 OpenCL 的 Command Queue，若沒有特別設定，原則上都是 in-order 處理 (相對於隨意順序，必須按照進入 queue 的順序執行)，因此 memory copy 也是一條指令，cudaMemcpy 屬於 blocking 函數，在設計上就不一定要加上 cudaDeviceSynchronize。

#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <string.h>
#include <cuda.h>
#define GPULOCAL 1024
#define CheckErr(status) { gpuAssert((status), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, const char *file, int line, int abort=true) {
    if (code != cudaSuccess) {
        fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
        if (abort) exit(code);
    }
}
#define MAXN 2048
static int N, M, BN;
static char str[MAXN][MAXN] = {};
static char hostMtx[2][MAXN*MAXN] = {};
__global__ void simulate(char IN[], char OUT[], int N, int BN) {
    int globalId = blockIdx.x * blockDim.x + threadIdx.x;
    int x = globalId/N + 1, y = globalId%N + 1;
#define G(x, y) IN[(x) * BN + (y)]
    char t = G(x, y);
    char adj = G(x-1, y-1) + G(x-1, y) + G(x-1, y+1)
            + G(x, y-1) + G(x, y+1) + G(x+1, y-1) + G(x+1, y) + G(x+1, y+1);
    OUT[x * BN + y] = (t == 0 && adj == 3) || (t == 1 && (adj == 2 || adj == 3));
#undef G
}
void runCuda() {
    int cudaDeviceCnt = 0;
    cudaGetDeviceCount(&cudaDeviceCnt);
    if (cudaDeviceCnt == 0) {
        printf("No supported GPU\n");
        return ;
    }
    
    char *cuMtx[2];
    int memSz = BN*BN*sizeof(char);
    int localSz = 1;
    for (int i = 1; i <= 1024; i++) {
        if (N*N%i == 0)
            localSz = i;
    }
    dim3 cuBlock(localSz);
    dim3 cuGrid(N*N/localSz);
    cudaMalloc((void **) &cuMtx[0], memSz);
    cudaMalloc((void **) &cuMtx[1], memSz);
    cudaMemcpy(cuMtx[0], hostMtx[0], memSz, cudaMemcpyHostToDevice);
    cudaMemcpy(cuMtx[1], hostMtx[1], memSz, cudaMemcpyHostToDevice);
    for (int i = 0; i < M; i++) {
        simulate<<<cuGrid, cuBlock>>>(cuMtx[i%2], cuMtx[(i+1)%2], N, BN);
        CheckErr(cudaGetLastError());
    }
    cudaDeviceSynchronize();	
    int f = M%2;
    cudaMemcpy(hostMtx[f], cuMtx[f], memSz, cudaMemcpyDeviceToHost);
    for (int i = 1; i <= N; i++) {
        for (int j = 1; j <= N; j++)
            hostMtx[f][i*BN + j] += '0';
        puts(hostMtx[f] + i*BN + 1);
    }
}
int main() {
    assert(scanf("%d %d", &N, &M) == 2);
    while (getchar() != '\n');
    for (int i = 1; i <= N; i++)
        assert(fgets(str[i]+1, MAXN, stdin) != NULL);
    BN = N+2;
    for (int i = 1; i <= N; i++) {
        for (int j = 1; j <= N; j++)
            hostMtx[0][i*BN + j] = str[i][j] - '0';
    }
    runCuda();
    return 0;
}

Morris' Blog