2016-06-28

批改娘 10100. Fast Matrix Multiplication (CUDA)

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

題目描述

計算兩個大小為 $N \times N$ 方陣 $A, \; B$ 相乘結果 $C = A \times B$。為了節省輸入輸出時間，採用亂數產生，可以參考下述程式碼，並改寫成 CUDA 的版本進行加速。

sequence.c

#include <stdio.h>
#include <stdint.h>
// #define DEBUG
#define UINT uint32_t
#define MAXN 1024
void multiply(int N, UINT A[][MAXN], UINT B[][MAXN], UINT C[][MAXN]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            UINT sum = 0;	// overflow, let it go.
            for (int k = 0; k < N; k++)
                sum += A[i][k] * B[k][j];
            C[i][j] = sum;
        }
    }
}
void rand_gen(UINT c, int N, UINT A[][MAXN]) {
    UINT x = 2, n = N*N;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            x = (x * x + c + i + j)%n;
            A[i][j] = x;
        }
    }
}
void print_matrix(int N, UINT A[][MAXN]) {
    for (int i = 0; i < N; i++) {
        fprintf(stderr, "[");
        for (int j = 0; j < N; j++)
            fprintf(stderr, " %u", A[i][j]);
        fprintf(stderr, " ]\n");
    }
}
UINT signature(int N, UINT A[][MAXN]) {
    UINT h = 0;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            h = (h + A[i][j]) * 2654435761LU;
    }
    return h;
}
UINT A[MAXN][MAXN], B[MAXN][MAXN], C[MAXN][MAXN];
int main() {
    int N;
    uint32_t S1, S2;
    scanf("%d %u %u", &N, &S1, &S2);
    rand_gen(S1, N, A);
    rand_gen(S2, N, B);
    multiply(N, A, B, C);
#ifdef DEBUG
    print_matrix(N, A);
    print_matrix(N, B);
    print_matrix(N, C);
#endif
    printf("%u\n", signature(N, C));
    return 0;
}

輸入格式

測資只有一組，包含三個整數 $N, S_1, S_2$，分別為方陣大小 $N \times N$，產生矩陣 $A$、$B$ 的亂數種子。

$64 \le N \le 1024$，保證 $N \mod 64 \equiv 0$
$0 \le S_1, \; S_2 < 2^{31}$

輸出格式

輸出一行雜湊值 $H$，可參考 sequence.c 的流程。

範例輸入

64 1 2

範例輸出

1	3376147904

編譯參數

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

Solution

與 OpenCL 版本相比，是否乾淨許多呢？完全不用管到底 Context 和 Kerenl Program 如何建造，但是這一些方便性都是因為預設導致的結果。

CUDA 預設會在 device 0 上運作，也就是 PCI-E 順位上的第一個顯卡，若要更動藉由函數 cudaSetDevice(deviceId); 完成，而 OpenCL 的 CommandQueue 則對應到 CUDA 的 Stream，在稍後的題目中，我會提供些許的範例使用 CUDA 這些函數，協助設計的計算流程可以加快。

#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <string.h>
#include <cuda.h>
#define MAXN 1024
#define GPULOCAL 64
#define UNLOOP 8
#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);
    }
}
uint32_t  hostA[MAXN*MAXN], hostB[MAXN*MAXN], hostC[MAXN*MAXN];
int N = MAXN;
__global__ void matrixMul(uint32_t A[], uint32_t B[], uint32_t C[], int N) {
    __shared__ uint32_t cbuf[MAXN+1];
    uint32_t rbuf[MAXN];
    int r = blockIdx.x * blockDim.x + threadIdx.x;
    int localId = threadIdx.x;
    int localSz = blockDim.x;
    for (int i = 0; i < N; i++)	
        rbuf[i] = A[r * N + i];
    for (int c = 0; c < N; c++) {
        for (int cr = localId; cr < N; cr += localSz)
            cbuf[cr] = B[cr * N + c];
        __syncthreads();
        uint32_t sum = 0;
        for (int k = 0; k+UNLOOP-1 < N; k += UNLOOP) {
            sum += rbuf[k+0] * cbuf[k+0];
            sum += rbuf[k+1] * cbuf[k+1];
            sum += rbuf[k+2] * cbuf[k+2];
            sum += rbuf[k+3] * cbuf[k+3];
            sum += rbuf[k+4] * cbuf[k+4];
            sum += rbuf[k+5] * cbuf[k+5];
            sum += rbuf[k+6] * cbuf[k+6];
            sum += rbuf[k+7] * cbuf[k+7];
        }
        C[r * N + c] = sum;
    }
}
void readIn() {
    uint32_t c1, c2;
    assert(scanf("%d %u %u", &N, &c1, &c2) == 3);
    uint32_t x = 2, n = N*N;
    x = 2;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            x = (x * x + c1 + i + j)&(n-1);
            hostA[i*N+j] = x;
        }
    }
    x = 2;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            x = (x * x + c2 + i + j)&(n-1);
            hostB[i*N+j] = x;
        }
    }
}
void writeOut() {
    uint32_t h = 0;
    uint32_t *Cend = hostC + N*N, *C = hostC;
    for (; C != Cend; C++)
        h = (h + *C) * 2654435761LU;
    printf("%u\n", h);
}
int main(int argc, char *argv[]) {
    readIn();
    uint32_t *cuMtxC, *cuMtxA, *cuMtxB;
    cudaMalloc((void **) &cuMtxC, N*N*sizeof(uint32_t));
    cudaMalloc((void **) &cuMtxA, N*N*sizeof(uint32_t));
    cudaMalloc((void **) &cuMtxB, N*N*sizeof(uint32_t));
    cudaMemcpy(cuMtxA, hostA, sizeof(uint32_t)*N*N, cudaMemcpyHostToDevice);	
    cudaMemcpy(cuMtxB, hostB, sizeof(uint32_t)*N*N, cudaMemcpyHostToDevice);
    CheckErr(cudaGetLastError());
    dim3 cuBlock(GPULOCAL);
    dim3 cuGrid(N/GPULOCAL);
    matrixMul<<<cuGrid, cuBlock>>>(cuMtxA, cuMtxB, cuMtxC, N);
    CheckErr(cudaGetLastError());
    cudaMemcpy(hostC, cuMtxC, sizeof(uint32_t)*N*N, cudaMemcpyDeviceToHost);
    CheckErr(cudaGetLastError());
    writeOut();	
    cudaFree(cuMtxC);
    return 0;
}

Morris' Blog