Day3: CUDA夏令营之矩阵乘法

已于 2022-07-06 17:28:53 修改
阅读量193 收藏
点赞数 1
分类专栏： CDUA基础 文章标签： 矩阵 线性代数
于 2022-07-06 17:14:32 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq906194732/article/details/125640426
版权

如何用cuda完成 矩阵相乘

一、 矩阵乘法 - 数学基础

  • 数学表达：
    在这里插入图片描述
  • 以3x3方阵乘法为例，搞清楚下标对应关系！
    在这里插入图片描述

注意：

  1. 矩阵A维度：m行 n列，矩阵B维度：j行 k列，二者能相乘的充要条件：n=j。即A矩阵的列数与B矩阵行数一致。
  2. 新矩阵维度 m行 k列。

二、矩阵乘法 - CPU实现

  • 本质：将矩阵坐标序列化到内存索引，分别申请3块内存用于计算。例：3x3矩阵序列化为 0~8。
  • 核心代码如下
。。。
void cpu_matrix_mult(int *h_a, int *h_b, int *h_result, int m, int n, int k) {
    // 遍历行
    /*
     A矩阵存在*h_a中，B矩阵存在*h_b中，新矩阵存在*h_result
     */
    for (int i = 0; i < m; ++i) 
    {
        // 遍历 列
        for (int j = 0; j < k; ++j) 
        {
            int tmp = 0.0;
            // 新矩阵中每个元素 需要进行 n次 乘加运算
            for (int h = 0; h < n; ++h) 
            {
                //索引计算：在A矩阵取第i行，B矩阵取第j列进行计算
                tmp += h_a[i * n + h] * h_b[h * k + j];  // 累加A矩阵第i行，B矩阵第j列乘积
            }
            h_result[i * k + j] = tmp;
        }
    }
}
。。。

三、矩阵乘法 - GPU实现

  • 难点：线程索引计算与线程分配。由于GPU的物理结构（Grid、Block、Thread），序列化矩阵坐标时，需要计算线程索引threadIdx，以保证矩阵相乘结果的正确性。
  • 内存申请
    int *h_a, *h_b, *h_c;
    cudaMallocHost((void **) &h_a, sizeof(int)*m*n); // 矩阵A
    cudaMallocHost((void **) &h_b, sizeof(int)*n*k); // 矩阵B
    cudaMallocHost((void **) &h_c, sizeof(int)*m*k); // A*B 结果

注意：此处函数后缀为Host，前缀为cuda，以 *h_a 为例

  1. Host代指CPU端，分配的是内存，大小为sizeof(int)*m *n，*h_a拿到内存地址
  2. cuda指，此函数与cuda有关，而非cuda控制CPU分配，Host是主体。
  • 显存申请
    int *d_a, *d_b, *d_c;
    cudaMalloc((void **) &d_a, sizeof(int)*m*n);
    cudaMalloc((void **) &d_b, sizeof(int)*n*k);
    cudaMalloc((void **) &d_c, sizeof(int)*m*k);
  • 数据拷贝 CPU -> GPU
    cudaMemcpy(d_a, h_a, sizeof(int)*m*n, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, sizeof(int)*n*k, cudaMemcpyHostToDevice);
  • 线程申请
    // 新矩阵 需要多少个线程
    // 向上取整，保证分配的 线程数 多余被拆分的 小问题
    unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
    unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;

    // 定义Grid维度，其中包含的 Grid 数量
    dim3 dimGrid(grid_cols, grid_rows);
    // 定义Block维度，其中包含的 Block 数量
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);

问题1：如何确定分配多少线程？
答：根据GPU计算时的最小子问题决定。最小子问题（可拆解为单一指令或者少指令）需要根据实际情况设定。例如矩阵乘法中，最小子问题为向量点乘，即计算新矩阵中每个位置的值。新矩阵维度为m * n，则至少需要分配 m * n个线程，用于 并行计算新矩阵每个位置的值。
问题2：如何确定GridDim与BlockDim的值？
答：由Block_size（每个Block含有的线程数目）与线程总数有关。Block_size设置有上限，一个Block中最多含有1024个线程。维度实际问题定。

// 仅使用 global memory 进行 矩阵乘法
__global__ void gpu_matrix_mult(int *a,int *b, int *c, int m, int n, int k)
// 矩阵A：m*n, 矩阵B：n*k,  *a: A,  *b: B, *c: C = A * B
{ 
    int row = blockIdx.y * blockDim.y + threadIdx.y; 
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int sum = 0;
    if( col < k && row < m) 
    {
        for(int i = 0; i < n; i++) 
        {
            sum += a[row * n + i] * b[i * k + col]; // 同 CPU
        }
        c[row * k + col] = sum;
    }
}

注意：此处只计算一行 * 一列，即，只定义单一指令操作，与CPU计算有很大不同。

  • 并行计算
    gpu_matrix_mult<<<dimGrid, dimBlock>>>(d_a, d_b, d_c, m,n,k);

此处为，调用线程进行并行计算，每个线程执行gpu_matrix_mult，相当于CPU中的两层for循环。

四、矩阵乘法 - GPU优化（使用Shared-memory）

  • Shared Memory:
    • 目前最快的多个线程可以交换数据的地方
    • 一个Block中的所有线程 共享一个Shared Memory
    • 逻辑上分为 32个bank，每个bank同一时间只能响应一个地址
    • 可以被设置为 16KB,32KB等

在这里插入图片描述

  • 如何利用Shared Memory来优化矩阵乘法呢？（三中使用的是Global memory）
// 使用 shared memory 进行 矩阵乘法 优化
__global__ void gpu_matrix_mult_shared(int *d_a, int *d_b, int *d_result, int n) 
{
    // shared memory 初始化
    __shared__ int tile_a[BLOCK_SIZE][BLOCK_SIZE];
    __shared__ int tile_b[BLOCK_SIZE][BLOCK_SIZE];

    // 行 = 前面有几行 * 每行元素数 + 这行的第几个
    int row = blockIdx.y * BLOCK_SIZE + threadIdx.y;
    // 列 = 前有几列 * 每列元素数 + 这列第几个
    int col = blockIdx.x * BLOCK_SIZE + threadIdx.x;
    int tmp = 0;
    int idx;

    // 将两个打矩阵相乘，拆分为 gridDim个 小矩阵对应相乘 再求和，以 sub 0表示：
    for (int sub = 0; sub < gridDim.x; ++sub) 
    {
        // 计算第 sub 部分矩阵的行索引
        idx = row * n + sub * BLOCK_SIZE + threadIdx.x;
        // 对 tile_a 进行赋值，即 global memory -> shared memory，多余部分置0
        tile_a[threadIdx.y][threadIdx.x] = row<n && (sub * BLOCK_SIZE + threadIdx.x)<n? d_a[idx]:0;
        // 计算第 sub 部分矩阵的行索引
        idx = (sub * BLOCK_SIZE + threadIdx.y) * n + col;
        // 对 tile_b 进行赋值，多余部分置0
        tile_b[threadIdx.y][threadIdx.x] = col<n && (sub * BLOCK_SIZE + threadIdx.y)<n? d_b[idx]:0;


        __syncthreads();  // 要同步，程序别退出
        // 计算 行、列 乘积，同时 对不同 sub 部分进行求和
        for (int k = 0; k < BLOCK_SIZE; ++k) 
        {
            tmp += tile_a[threadIdx.y][k] * tile_b[k][threadIdx.x];
        }
        __syncthreads();  // 要同步，程序别退出
    }
    // 那么d_result中序号row * n + col的元素值 就为 tmp
    if(row < n && col < n)
    {
        d_result[row * n + col] = tmp;
    }
}

五、完整代码理解及测试结果

  • 使用 nvida-nona进行测试，结果如下：
    在这里插入图片描述

结论：使用GPU进行矩阵相乘计算，真快！

  • 源代码如下：
#define CHECK(call)                                   \
do                                                    \
{                                                     \
    const cudaError_t error_code = call;              \
    if (error_code != cudaSuccess)                    \
    {                                                 \
        printf("CUDA Error:\n");                      \
        printf("    File:       %s\n", __FILE__);     \
        printf("    Line:       %d\n", __LINE__);     \
        printf("    Error code: %d\n", error_code);   \
        printf("    Error text: %s\n",                \
            cudaGetErrorString(error_code));          \
        exit(1);                                      \
    }                                                 \
} while (0)

#include <stdio.h>
#include <math.h>
#include "error.cuh"

#define BLOCK_SIZE 16

// 仅使用 global memory 进行 矩阵乘法
__global__ void gpu_matrix_mult(int *a,int *b, int *c, int m, int n, int k)
// 矩阵A：m*n, 矩阵B：n*k,  *a: A,  *b: B, *c: C = A * B
{ 
    int row = blockIdx.y * blockDim.y + threadIdx.y; 
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int sum = 0;
    if( col < k && row < m) 
    {
        for(int i = 0; i < n; i++) 
        {
            sum += a[row * n + i] * b[i * k + col];
        }
        c[row * k + col] = sum;
    }
} 

// 使用 shared memory 进行 矩阵乘法 优化
__global__ void gpu_matrix_mult_shared(int *d_a, int *d_b, int *d_result, int n) 
{
    // shared memory 初始化
    __shared__ int tile_a[BLOCK_SIZE][BLOCK_SIZE];
    __shared__ int tile_b[BLOCK_SIZE][BLOCK_SIZE];

    // 行 = 前面有几行 * 每行元素数 + 这行的第几个
    int row = blockIdx.y * BLOCK_SIZE + threadIdx.y;
    // 列 = 前有几列 * 每列元素数 + 这列第几个
    int col = blockIdx.x * BLOCK_SIZE + threadIdx.x;
    int tmp = 0;
    int idx;

    // 将两个打矩阵相乘，拆分为 gridDim个 小矩阵对应相乘 再求和，以 sub 0表示：
    for (int sub = 0; sub < gridDim.x; ++sub) 
    {
        // 计算第 sub 部分矩阵的行索引
        idx = row * n + sub * BLOCK_SIZE + threadIdx.x;
        // 对 tile_a 进行赋值，即 global memory -> shared memory，多余部分置0
        tile_a[threadIdx.y][threadIdx.x] = row<n && (sub * BLOCK_SIZE + threadIdx.x)<n? d_a[idx]:0;
        // 计算第 sub 部分矩阵的行索引
        idx = (sub * BLOCK_SIZE + threadIdx.y) * n + col;
        // 对 tile_b 进行赋值，多余部分置0
        tile_b[threadIdx.y][threadIdx.x] = col<n && (sub * BLOCK_SIZE + threadIdx.y)<n? d_b[idx]:0;


        __syncthreads();  // 要同步，程序别退出
        // 计算 行、列 乘积，同时 对不同 sub 部分进行求和
        for (int k = 0; k < BLOCK_SIZE; ++k) 
        {
            tmp += tile_a[threadIdx.y][k] * tile_b[k][threadIdx.x];
        }
        __syncthreads();  // 要同步，程序别退出
    }
    // 那么d_result中序号row * n + col的元素值 就为 tmp
    if(row < n && col < n)
    {
        d_result[row * n + col] = tmp;
    }
}

// 使用 cpu 进行 矩阵乘法
void cpu_matrix_mult(int *h_a, int *h_b, int *h_result, int m, int n, int k) {
    // 遍历行
    /*
     A矩阵存在*h_a中，B矩阵存在*h_b中，新矩阵存在*h_result
     */
    for (int i = 0; i < m; ++i) 
    {
        // 遍历 列
        for (int j = 0; j < k; ++j) 
        {
            int tmp = 0.0;
            // 新矩阵中每个元素 需要进行 n次 乘加运算
            for (int h = 0; h < n; ++h) 
            {
                //索引计算：在A矩阵取第i行，B矩阵取第j列进行计算
                tmp += h_a[i * n + h] * h_b[h * k + j];  // 累加A矩阵第i行，B矩阵第j列乘积
            }
            h_result[i * k + j] = tmp;
        }
    }
}

int main(int argc, char const *argv[])
{
    // 矩阵维度 m*n n*k -> m*k
    int m=1000;
    int n=1000;
    int k=1000;

    // 分配 内存 并检查
    int *h_a, *h_b, *h_c, *h_cc, *h_cs;
    CHECK(cudaMallocHost((void **) &h_a, sizeof(int)*m*n));
    CHECK(cudaMallocHost((void **) &h_b, sizeof(int)*n*k));
    CHECK(cudaMallocHost((void **) &h_c, sizeof(int)*m*k));
    CHECK(cudaMallocHost((void **) &h_cc, sizeof(int)*m*k));
    CHECK(cudaMallocHost((void **) &h_cs, sizeof(int)*m*k));

    // 创建 event 用于统计时间
    cudaEvent_t start, stop,stop_share;
    CHECK(cudaEventCreate(&start));
    CHECK(cudaEventCreate(&stop));
    CHECK(cudaEventCreate(&stop_share));


    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            h_a[i * n + j] = 1;
        }
    }

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < k; ++j) {
            h_b[i * k + j] = 0;
        }
    }

    int *d_a, *d_b, *d_c, *d_c_share;
    // 分配显存 并检查
    CHECK(cudaMalloc((void **) &d_a, sizeof(int)*m*n));
    CHECK(cudaMalloc((void **) &d_b, sizeof(int)*n*k));
    CHECK(cudaMalloc((void **) &d_c, sizeof(int)*m*k));
    CHECK(cudaMalloc((void **) &d_c_share, sizeof(int)*m*k));

    CHECK(cudaEventRecord(start));
    // copy matrix A and B from host to device memory
    CHECK(cudaMemcpy(d_a, h_a, sizeof(int)*m*n, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(d_b, h_b, sizeof(int)*n*k, cudaMemcpyHostToDevice));

    // 新矩阵 需要多少个线程
    // 向上取整，保证分配的 线程数 多余被拆分的 小问题
    unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
    unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;

    // 定义Grid维度，其中包含的 Grid 数量
    dim3 dimGrid(grid_cols, grid_rows);
    // 定义Block维度，其中包含的 Block 数量
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
    
    // 实际执行 use global memor
    gpu_matrix_mult<<<dimGrid, dimBlock>>>(d_a, d_b, d_c, m,n,k);    

    // 将新 矩阵，GPU -> CPU
    CHECK(cudaMemcpy(h_c, d_c, (sizeof(int)*m*k), cudaMemcpyDeviceToHost));
    //cudaThreadSynchronize();
    // 结束 stop event
    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));


    // 实际执行 use shared memor
    gpu_matrix_mult_shared<<<dimGrid, dimBlock>>>(d_a, d_b, d_c_share, n);
    CHECK(cudaMemcpy(h_cs, d_c_share, (sizeof(int)*m*k), cudaMemcpyDeviceToHost));

    // 结束 stop stop_share
    CHECK(cudaEventRecord(stop_share));
    CHECK(cudaEventSynchronize(stop_share));
    
    float elapsed_time, elapsed_time_share;
    CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));
    CHECK(cudaEventElapsedTime(&elapsed_time_share, stop, stop_share));
    printf("Time_global = %g ms.\n", elapsed_time);
    printf("Time_share = %g ms.\n", elapsed_time_share);

    CHECK(cudaEventDestroy(start));
    CHECK(cudaEventDestroy(stop));    

    // 实际执行 CPU
    cpu_matrix_mult(h_a, h_b, h_c, m, n, k);

    int ok = 1;
    for (int i = 0; i < m; ++i)
    { 
        for (int j = 0; j < k; ++j)
        {
            if(fabs(h_cs[i*k + j] - 0)>(1.0e-10))
            {
                printf("hcs: %d hc: %d  ",h_cs[i*k + j], h_c[i*k + j]);
                ok = 0;
            }
        }
    }

    if(ok)
    {
        printf("Pass!!!\n");
    }
    else
    {
        printf("Error!!!\n");
    }
    
    // free memory
    CHECK(cudaFree(d_a));
    CHECK(cudaFree(d_b));
    CHECK(cudaFree(d_c));
    CHECK(cudaFreeHost(h_a));
    CHECK(cudaFreeHost(h_b));
    CHECK(cudaFreeHost(h_c));
    return 0;
}
passioq_
关注
  • 1
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值