cuda是如何计算的

矩阵乘法是啥，在这里就不多说了。

当我们的block dim是11的时候MatmulKernel函数每次会计算一个位置的结果。
当block dim是22的时候MatmulKernel函数每次会计算4个位置的结果。
无论block dim多大，每个线程都只计算一个位置的结果。

#include "matmul_basic.h"
#include "cuda_runtime.h"
#include "cuda.h"
#include "stdio.h"

/* matmul的函数实现*/
__global__ void MatmulKernel(float *M_device, float *N_device, float *P_device, int width){
    /* 
        我们设定每一个thread负责P中的一个坐标的matmul
        所以一共有width * width个thread并行处理P的计算
    */
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    int x = blockIdx.x * blockDim.x + threadIdx.x;

    float P_element = 0;

    /* 对于每一个P的元素，我们只需要循环遍历width次M和N中的元素就可以了*/
    for (int k = 0; k < width; k ++){
        float M_element = M_device[y * width + k];
        float N_element = N_device[k * width + x];
        P_element += M_element * N_element;
    }

    P_device[y * width + x] = P_element;
}

void MatmulOnDevice(float *M_host, float *N_host, float* P_host, int width, int blockSize){
    /* 设置矩阵大小 */
    int size = width * width * sizeof(float);

    /* 分配M, N在GPU上的空间*/
    float *M_device;
    float *N_device;

    cudaMalloc(&M_device, size);
    cudaMalloc(&N_device, size);

    /* 分配M, N拷贝到GPU上*/
    cudaMemcpy(M_device, M_host, size, cudaMemcpyHostToDevice);
    cudaMemcpy(N_device, N_host, size, cudaMemcpyHostToDevice);

    /* 分配P在GPU上的空间*/
    float *P_device;
    cudaMalloc(&P_device, size);

    /* 调用kernel来进行matmul计算：将一个矩阵切分成多个blockSize * blockSize的大小 */
    dim3 dimBlock(blockSize, blockSize);
    dim3 dimGrid(width / blockSize, width / blockSize);  //有width*width个block
    MatmulKernel <<<dimGrid, dimBlock>>> (M_device, N_device, P_device, width);

    /* 将结果从device拷贝回host*/
    cudaMemcpy(P_host, P_device, size, cudaMemcpyDeviceToHost);
    cudaDeviceSynchronize();

    /* Free */
    cudaFree(P_device);
    cudaFree(N_device);
    cudaFree(M_device);
}

完整代码

// matmul_basic.h
#pragma once

void MatmulOnDevice(float *M_host, float *N_host, float* P_host, int width, int blockSize);

//main.cpp
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include "timer.hpp"
#include "matmul_basic.h"
#include "util.h"

void MatmulOnHost(float *M, float *N, float *P, int width){
    for (int i = 0; i < width; i ++)
        for (int j = 0; j < width; j ++){
            float sum = 0;
            for (int k = 0; k < width; k++){
                float a = M[i * width + k];
                float b = N[k * width + j];
                sum += a * b;
            }
            P[i * width + j] = sum;
        }
}


int main(){
    int seed;
    Timer timer;
    int width  = 256;
    int min  = 0;
    int max = 5;
    int size_matrix = width * width;
    int block_width = 1;

    float* h_matM = (float*)malloc(size_matrix * sizeof(float));
    float* h_matN = (float*)malloc(size_matrix * sizeof(float));
    float* h_matP = (float*)malloc(size_matrix * sizeof(float));
    float* d_matP = (float*)malloc(size_matrix * sizeof(float));
    
    //生成2个随机Matrix
    seed = 1;
    initMatrix(h_matM, size_matrix, min, max, seed);
    seed += 1;
    initMatrix(h_matN, size_matrix, min, max, seed);
    
    /* CPU */
    timer.start();
    MatmulOnHost(h_matM, h_matN, h_matP, width);
    // print_Matrix(h_matP, size_matrix);
    timer.stop();
    timer.duration<Timer::ms>("matmul in cpu");

    /* GPU warmup */
    timer.start();
    MatmulOnDevice(h_matM, h_matN, d_matP, width, block_width);
    timer.stop();
    timer.duration<Timer::ms>("matmul in gpu(warmup)");

    timer.start();
    MatmulOnDevice(h_matM, h_matN, d_matP, width, block_width);
    timer.stop();
    timer.duration<Timer::ms>("matmul in gpu(bs = 1)");
     compareMat(h_matP, d_matP, size_matrix);

    /* GPU general implementation, bs = 16*/
    block_width = 16;
    timer.start();
    MatmulOnDevice(h_matM, h_matN, d_matP, width, block_width);
    timer.stop();
    timer.duration<Timer::ms>("matmul in gpu(bs = 16)");
    compareMat(h_matP, d_matP, size_matrix);

    // /* GPU general implementation, bs = 32*/
    block_width = 32;
    timer.start();
    MatmulOnDevice(h_matM, h_matN, d_matP, width, block_width);
    timer.stop();
    timer.duration<Timer::ms>("matmul in gpu(bs = 32)");
    compareMat(h_matP, d_matP, size_matrix);

    /* 
     * 注意，这里blockSize=64导致一个block中的thread数量超过了1024，最终使得kernel无法启动
     * 这个错误属于参数设定的错误。类似的错误比如说还有设置过大的shared_memory
    */
    timer.start();
    block_width = 64;
    MatmulOnDevice(h_matM, h_matN, d_matP, width, block_width);
    timer.stop();
    // std::sprintf(str, "matmul in gpu(general)<<<%d, %d>>>", width / blockSize, blockSize);
    timer.duration<Timer::ms>("matmul in gpu(bs = 64)");
    compareMat(h_matP, d_matP, size_matrix);
    // return 0;
}

//util.cpp
#include "util.h"
#include <stdlib.h>
#include <cstdio>


void initMatrix(float* data, int size, int min, int max, int seed) {
    srand(seed);
    for (int i = 0; i < size; i ++) {
        data[i] = float(rand()) * float(max - min) / RAND_MAX;
    }
}

void printMat(float* data, int size) {
    for (int i = 0; i < size; i ++) {
        printf("%.8lf", data[i]);
        if (i != size - 1) {
            printf(", ");
        } else {
            printf("\n");
        }
    }
}

void compareMat(float* h_data, float* d_data, int size) {
    //浮点数运算时CPU和GPU之间的计算结果是有误差的
    double precision = 1.0E-3;

    for (int i = 0; i < size; i ++) {
        if (abs(h_data[i] - d_data[i]) > precision) {
            int y = i / size;
            int x = i % size;
            printf("Matmul result is different\n");
            printf("cpu: %.8lf, gpu: %.8lf, cord:[%d, %d]\n", h_data[i], d_data[i], x, y);
            break;
        }
    }
}

//util.h
#pragma once

void initMatrix(float* data, int size, int low, int high, int seed);

void printMat(float* data, int size);

void compareMat(float* h_data, float* d_data, int size);

//timer.hpp
#include <chrono>
#include <cstdio>
#include <ratio>
#include <string>
#include <iostream>


class Timer {
public:
    using s  = std::ratio<1, 1>;
    using ms = std::ratio<1, 1000>;
    using us = std::ratio<1, 1000000>;
    using ns = std::ratio<1, 1000000000>;

public:
    Timer(){};

public:
    void start() {mStart = std::chrono::high_resolution_clock::now();}
    void stop()  {mStop  = std::chrono::high_resolution_clock::now();}

    template <typename T>
    void duration(std::string msg);

private:
    std::chrono::time_point<std::chrono::high_resolution_clock> mStart;
    std::chrono::time_point<std::chrono::high_resolution_clock> mStop;
};

template <typename T>
void Timer::duration(std::string msg){
    std::string str;
    char fMsg[100];
    std::sprintf(fMsg, "%-30s", msg.c_str());

    if(std::is_same<T, s>::value) { str = " s"; }
    else if(std::is_same<T, ms>::value) { str = " ms"; }
    else if(std::is_same<T, us>::value) { str = " us"; }
    else if(std::is_same<T, ns>::value) { str = " ns"; }

    std::chrono::duration<double, T> time = mStop - mStart;
    std::cout << fMsg << " uses " << time.count() << str << std::endl;
}

#CMakeLists.txt
cmake_minimum_required(VERSION 2.8) # cmake最低版本
project(manager_project CUDA CXX) #项目名称
set(CMAKE_CXX_STANDARD 11) #设置C++编译版本
set(CMAKE_CUDA_STANDARD 11)
set(CMAKE_BUILD_TYPE "Debug") # 默认是Release模式，设置为Debug才能调试
set(EXECUTABLE_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/bin) #设置可执行文件生产的路径
file(GLOB SOURCES ${PROJECT_SOURCE_DIR}/*.cpp ${PROJECT_SOURCE_DIR}/*.cu )
add_executable(demo ${SOURCES} ) #生成可执行文件demo