【参加CUDA线上训练营】储存单元矩阵乘法

知识点

上节课问题部分答疑

1. 一个block只有1024个线程。
2. 最佳参数需要自己尝试
3. 每个并行程序输出乱序，是因为输出buff导致。
4. 注意判断矩阵数据越界。

存储单元

与CPU（host）交互的内存Global、Constant、Texture Memory
此外线程可与其他的shared memory、 register、local memory等内存单元读写。

内存申请操作

GPU与CPU数据异步运行，且需要提前在GPU上申请数据空间，主要功能函数为cudaMemcpy()、cudaFree()。
gpu和cpu异步运行，需要同步函数

矩阵相乘

cpu：

CPU上计算矩阵乘法需要重复计算乘法并累加，计算繁杂，花费较大

 void cpu_matrix_mult(int *h_a, int *h_b, int *h_result, int m, int n, int k) 
    {
        for (int i = 0; i < m; ++i) 
        {
            for (int j = 0; j < k; ++j) 
            {
                int tmp = 0;
                for (int h = 0; h < n; ++h) 
                {
                    tmp += h_a[i * n + h] * h_b[h * k + j];
                }
                h_result[i * k + j] = tmp;
            }
        }
    }

如果矩阵a，b维度1000，就要作1000^3次tmp += h_a[i * n + h] * h_b[h * k + j];

GPU

利用矩阵索引替代for循环中的i，j
矩阵乘法=每个元素乘积累加的运算。每个无关元素运算互不影响
GPU提高速度：利用空间换取时间。多个核心替代运算步骤

实际上，二维矩阵数据储存在一维的连续空间上，所以可以根据x，y坐标找到二维矩阵中的元素。
index索引计算示例如下图。threadidx(2,0)为网格单元中所示单元数字, blockidx(1,1)为所示单元所在的block坐标, blockDim(4,3)为每个block的大小。

实验实例

矩阵乘法示例
![
对于 sum+=M[row*n_size+i]*N[i*k_size+col];
M是行乘，所有元素内存位置相邻，故是

行数*宽度（起始位置）+i（当行迭代索引）

N是列乘，每个元素内存位置都相差一个行元素的位置，故是

i*列数（起始位置）+col（高度），内存中离散位置

#include <stdio.h>
#include <cuda_runtime.h>
#include "device_launch_parameters.h"
#include <math.h>
#define BLOCK_SIZE 16
__global__ void gpu_matrix_mult(int *M,int *N, int *P,int m_size,int n_size, int k_size){
    //计算当前线程坐标
    int row=blockIdx.x*blockDim.x+threadIdx.x;
    int col=blockIdx.y*blockDim.y+threadIdx.y;

    //乘积运算
    int sum=0;
    if(row<m_size&&col<k_size){
        for(int i=0;i<n_size;i++){
            sum+=M[row*n_size+i]*N[i*k_size+col];
        }
        P[row*k_size+col]=sum;
    }
}

int main(int argc, char const *argv[])
{
    int m=100;
    int n=100;
    int k=100;
    // 申请计算内存
    int *h_a, *h_b, *h_c, *h_cc;
    cudaMallocHost((void **) &h_a, sizeof(int)*m*n);
    cudaMallocHost((void **) &h_b, sizeof(int)*n*k);
    cudaMallocHost((void **) &h_c, sizeof(int)*m*k);
    cudaMallocHost((void **) &h_cc, sizeof(int)*m*k);

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

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < k; ++j) {
            h_b[i * k + j] = rand() % 1024;
        }
    }
    //申请结果内存
    int *d_a, *d_b, *d_c;
    cudaMalloc((void**)&d_a,sizeof(int)*m*n);
    cudaMalloc((void**)&d_b,sizeof(int)*k*n);
    cudaMalloc((void**)&d_c,sizeof(int)*m*k);

    //转移数据
    cudaMemcpy(d_a,h_a,sizeof(int)*m*n,cudaMemcpyHostToDevice);
    cudaMemcpy(d_b,h_b,sizeof(int)*k*n,cudaMemcpyHostToDevice);

    //计算
    unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
    unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;
    dim3 dimGrid(grid_cols, grid_rows);
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);

    gpu_matrix_mult<<<dimGrid,dimBlock>>>(d_a,d_b,d_c,m,n,k);

    cudaMemcpy(h_c,d_c,sizeof(int)*m*k,cudaMemcpyDeviceToHost);

        // free memory
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    cudaFreeHost(h_a);
    cudaFreeHost(h_b);
    cudaFreeHost(h_c);
    cudaFreeHost(h_cc);
    return 0;
}

利用nvprof查看
size=100
size=1000

总结

1. 两个指针是指针的指针，cuda内存用指针指向数据，如要修改指针故要用两个星号。
2. gpu和cpu异步运行，需要同步函数
3. blockid要根据grid的分块判断，上图就把grid分为一个块，所以整个idx=(0,0)，Dim=(12,6),threadidx
   =(6,3),所以整个结果为(0,0)*(12,6)+(6,3)=(6,3)
   
4. 通过工具分析，其实增加size=1000后计算所用时间变化并不大，说明了GPU并行计算的高校性。
5. GPU大部分时间花在计算上，Host到Device的时间相比device到Host花销更大。