【CUDA编程3】简单的矩阵乘法

CUDA的优势是并行计算，计算机中最常见的需要并行计算的地方就是矩阵运算。图像图像处理、数字信号处理、神经网络算法等等都包含大量的矩阵运算。这也是CUDA被广泛使用的原因。
了解完CUDA程序的基本构成、CUDA程序的资源分配之后，可以实战进行一个矩阵乘法。
Demo如下：

#include<stdio.h>
#include<cuda.h>
#include<cuda_runtime.h>

#define BLOCK_NUM 4  
#define THREAD_NUM 4
#define R_SIZE BLOCK_NUM * THREAD_NUM
#define M_SIZE R_SIZE * R_SIZE

__global__ void mat_mul(int *mat1, int *mat2, int *result) {
    const int bid = blockIdx.x;
    const int tid = threadIdx.x;
    const int row = bid * THREAD_NUM + tid;
    for (int c = 0; c < R_SIZE; c++) {
        for (int n = 0; n < R_SIZE; n++) {
            result[row*R_SIZE+c] += mat1[row*R_SIZE+n] * mat2[n*R_SIZE+c];
        }
    }
}

int main(int argc, char *argv[]) {
    int *mat1, *mat2, *result;
    int *g_mat1, *g_mat2, *g_mat_result;
    
    // 1-dim NxN vector to represent 2-dim (N, N) matrix
    mat1 = (int*) malloc(M_SIZE * sizeof(int));
    mat2 = (int*) malloc(M_SIZE * sizeof(int));
    result = (int*) malloc(M_SIZE * sizeof(int));
    printf("M_SIZE:%d\n", M_SIZE);
    // init matrices
    for (int i = 0; i < M_SIZE; i++) {
        mat1[i] = rand() % 10;
        mat2[i] = rand() % 10;
        result[i] = 0;
    }
    cudaMalloc((void **)&g_mat1, sizeof(int) * M_SIZE);
    cudaMalloc((void **)&g_mat2, sizeof(int) * M_SIZE);
    cudaMalloc((void **)&g_mat_result, sizeof(int) * M_SIZE);
    cudaMemcpy(g_mat1, mat1, sizeof(int) * M_SIZE, cudaMemcpyHostToDevice);
    cudaMemcpy(g_mat2, mat2, sizeof(int) * M_SIZE, cudaMemcpyHostToDevice);
    mat_mul<<<BLOCK_NUM, THREAD_NUM>>>(g_mat1, g_mat2, g_mat_result);
    cudaMemcpy(result, g_mat_result, sizeof(int) * M_SIZE, cudaMemcpyDeviceToHost);
    printf("res[0]:%d\n", result[0]);
}

保存为matmul.cu，然后编译运行：

nvcc -o matmul matmul.cu
./matmul

这里用了4个block，每个block含有4个thread。处理2个16x16大小的矩阵乘法。相当于用16个线程并行计算这个矩阵乘法。

首先把16x16大小的两个矩阵都用256x1的向量来表示。通过用向量索引的二维性表示来实现矩阵的表示。