cuda开发基础

cuda开发基础

1 CUDA简介

CUDA(Compute Unified Device Architecture)的中文全称为计算同样设备架构。在高性能计算领域，CUDA是一个重要的工具。

2 GPU架构特点

首先，谈一谈串行计算和并行计算。高性能计算领域的关键是领域多核处理器进行并行计算。

当求解一个计算机程序任务时，很自然的想法是将任务分解成一系列小任务，然后把这些任务一一完成。在串行计算时，只能让处理器每次处理一个计算任务，处理完一个计算任务后再计算下一个任务，直到所有小任务都完成后，这个大任务就完成了。如下图所示，即为用串行编程思想求解问题的步骤。

 然后，串行计算的缺点非常明显，当我们拥有多核处理器，可以利用多核处理器同时处理多个任务时，且这些小任务并没有关联关系，那么就没必要使用串行编程的求解思路了。为了进一步加快大任务的计算速度，我们可以把一些独立的模块分配到不同处理器上进行同时计算（这就是并行），最后再将这些结果进行整合，完成一次任务计算。下图就是将一个大的计算任务分解为小任务，然后将独立的小任务分配到不同处理器进行并行计算，最后通过串行程序把结果汇总完成这次总的计算任务。

 因此，一个程序可不可以进行并行计算，关键在于程序可以拆分为那几个执行模块，这些执行模块哪些是独立的，哪些是强依赖强耦合的，独立的模块则可以设计为并行计算，并充分利用多核处理器优势进一步加速我们的计算任务，强耦合模块则使用串行编程。利用串行+并行的编程思想完成一次高性能计算。

再谈CPU和GPU的区别和各自特点。在并行、串行计算时多次谈到“多核”的概念。那么先从“核”的角度开始这个话题。首先CPU是专为顺序串行处理而优化的几个核心组成。而GPU则由数以千计的更小、更高效的核心组成。这些核心专门为同时处理多任务而设计，可以高效地处理并行任务。也就是CPU虽然每个核心的逻辑处理能力比较强，但本身核心较少，不适用于处理矩阵运算。但是GPU虽然每个核心逻辑处理能力很弱，但是由于计算核心很多，适合于大矩阵运算。这使得在处理多个计算任务时，在并行计算的支持上做得更高效。比如在计算机视觉处理领域，神经网络中的卷积运算。

GPU和CPU硬件的特点决定了他们的应用场景，CPU是计算机运算和控制的核心，GPU主要用于图像图像处理。图像在计算机呈现的形式就是矩阵，我们对图像的处理其实就是操作各种矩阵进行计算，而很多矩阵的运算其实可以做并行化，这使得图像处理可以做得很快，因此GPU在图像领域也就有了大展拳脚的机会。下图表示一个多GPU计算机硬件系统，可以看出，一个GPU内存就有很多个SP和各类内存。这些硬件都市GPU进行高效并行计算的基础。

 现在再从数据处理的角度来对比CPU和GPU的特点。CPU需要很强的通用性来处理各种不同的数据类型，比如整型、浮点型等，同时它又必须擅长处理逻辑判断所导致的大量分支和中断处理。而GPU面对的则是类型高度统一，相互无依赖的大规模数据和不需要被打断的纯净的计算环境，GPU有非常多核心（费米架构就有512核），虽然其核心能力远没有CPU的核心强，但是胜在多。再处理简单计算任务时呈现“人多力量大”的优势，这就是并行计算的魅力。

所以二者特点分别是：

• CPU：擅长流程控制和逻辑处理，不规则数据结构，不可预测存储结构，单线程程序，分支密集型算法。
• GPU：擅长数据并行处理，规则数据结构，可预测存储模式。

 在现在计算机体系架构中，要完成CUDA并行计算需要CPU和GPU协同配合才能完成一次高性能计算任务。

一般而言，并行部分在GPU上运行，串行部分在CPU上运行，这就是异构计算。异构计算就是不同体系结构的处理器互相协作完成一次计算任务。CPU负责总体程序流程，而GPU负责具体的计算任务，当GPU个线程完成计算任务后，我们就将GPU计算得到的结果拷贝到CPU端，完成一次计算任务。

 因此，应用程序利用GPU实现加速的总体分工：密集计算代码由GPU负责执行，剩余的串行代码由CPU负责执行。

3 CUDA线程模型

线程是程序执行的最基本单元，CUDA并行计算就是通过成千上万个线程的并行执行来实现的。下图说明了GPU的线程模型。

 

 CUDA的线程模型从小往大依次是：

Thread：线程，并行的基本单位。

Thread Block：线程块，互相合作的线程组，线程块有如下几个特点：

• 允许彼此同步
• 可以通过共享内存快速交换数据
• 以1维、2维、3维组织

Grid：一组线程块

• 以1维、2维组织
• 共享全局内存

Kernel：在GPU上执行的核心程序，这个kernel函数是运行在某个Grid上的。

• One kernel<->One Grid

每一个block和thread都有自己的ID，我们通过相应的索引找到相应的线程和线程块。

• threadidx、blockidx
• Block ID：1D or 2D
• Thread ID：1D、2D or 3D

理解kernel，必须对kernel的线程层次结构有一个清晰的认识。首先GPU上很多并行化的轻量级线程。kernel在device上执行时实际上是启动很多线程，一个kernel所启动的所有线程称为一个网格（Grid），同一个网格上的线程共享相同的全局内存空间，Grid是线程结构的第一个层次，而网格又可以分为很多线程块（Block），一个线程块中包含很多线程，这是线程结构的第二个层次。线程两层组织结构如上图所示，这个一个Grid和Block均为2-Dim的线程组织。Grid和Block都是定义为Dim3类型的变量，dim3可以看成是一个包含三个无符号整数(x,y,z)成员的结构体变量，在定义时，缺省值初始化为1。因此grid和block可以灵活地定义为1-dim，2-dim，3-dim结构，kernel调用时也必须通过配置<<<grid,block>>>来指定kernel所使用的网格维度和线程块维度。举个例子，我们以上图为例，分析怎么通过<<<grid,block>>>这种标记方式索引到我们想要的那个线程。CUDA的这种<<<grid,block>>>其实就是一个多级索引的方法，第一级索引是(grid.xidx,grid,.yidy)，对应上面例子就是(1,1)，通过它我们就能从找到这个线程块的位置，然后启动二级索引(block.xIdx, block.yIdx, block.zIdx)来定位到指定的线程。这就是CUDA的线程组织结构。

• SP：最基本的处理单元，streaming processor，也称为CUDA core。具体的指令和任务都是在SP上处理。GPU进行并行计算，也就是很多个SP同时做处理。 
• SM：多个SP加上其他一些资源组成一个streaming mutiprocessor。也叫GPU大核，其他资源如：warp scheduler，register，shared memory等。SM可以看作GPU的心脏，register和shared memory是SM的资源。CUDA将这些资源分配给所有驻留在SM中的threads。因此，这些有限的资源就使每个SM中的active warps有非常严格的限制，也就限制了并发能力。

需要指出，每个SM包含的SP数量依据GPU架构而不同，Fermi架构GF100是32个，GF10X是48个，Kepler架构都是192个，Maxwell都是128个。

简而言之，SP是线程执行的硬件单位，SM中包含多个SP，一个GPU可以有多个SM（比如16个），最终一个GPU可以能包含有上千个SP。这么多核心“是同时运行”，速度可想而知。这个引号只是想表明实际上，软件逻辑上是多有SP是并行的，但是物理上并不是所有SP都同时执行计算（比如我们只有8个SM却有1024个线程块需要调度处理），因为有些会处于挂起，就绪等其他状态，这有关GPU的线程调度。

下面这个图将从硬件和软件角度解释CUDA的线程模型。

•  每个线程由每个线程处理器（SP）执行
• 线程块由多核处理器（SM）执行
• 一个kernel其实由一个grid来执行，一个kernel一次只能一个GPU上执行

block是软件概念，一个block只会由一个sm调度。程序员在开发时，通过设定block属性，告诉GPU硬件，我有多少个线程，线程怎么组织。而具体怎么调度有sm的warps scheduler负责，block一旦被分配好sm，该block就会一直驻留在该SM中，直到执行结束。一个SM可以同时拥有多个blocks，但需要序列执行。下图显示了GPU内部硬件架构：

 4 CUDA内存模型

CUDA的内存模型分为以下几个层次：

• 每个线程独有自己的registers（寄存器）
• 每个线程独有自己的local memory（局部内存）
• 每个线程块独有自己的shared memory（共享内存），所有线程块内的所有线程共享这段内存资源
• 每个grid独有自己的global memory（全局内存），不同线程块的线程都可以使用
• 每个grid都有自己的constant memory（常量内存）和texture memory（纹理内存），不同线程块的线程都可以使用

线程访问这几类存储器的速度：registers>local memory>shared memory>global memory

下图表示的就是这些内存在计算机架构中的所在层次。

 5 CUDA编程模型

 上面讲了这么多硬件相关的知识点，现在来使用CUDA编写程序。

先捋一捋场景的CUDA术语：

 5.1 第一个编程要点

我们怎么编写一个能在GPU上跑的程序或函数呢？

通过关键字就可以表示某个程序在GPU上运行还是在CPU上运行。如下表所示，比如用__global__定义一个kernel函数，就是在CPU上调用，GPU上执行，注意__global__函数的返回值必须设置为void。

 5.2 第二个编程要点

CPU和GPU间的数据如何传输？

首先介绍在GPU内存分配和内存回收函数接口：

• cudaMalloc()：在设备端分配global memory
• cudaFree()：释放存储空间

CPU端数据和GPU端数据传输的函数接口时一样的，他们通过传递函数实参来表示传输方向：

• cudaMemcpy(void *dst,void *src,size_t nbytes,enum cudaMemcpyKind direction)
  • cudaMemcpyHostToDevice（CPU到GPU）
  • cudaMemcpyDeviceToHost（GPU到CPU）
  • cudaMemcpyDeviceToDevice（GPU到GPU）

5.3 第三个编程要点

 怎么用代码表示线程组织模型？

我们可以用dim3类来表示网格和线程块的组织方式，网格grid可以表示为一维和二维格式，线程块block可以表示维一维、二维和三维的数据格式。

dim3 DimGrid(100, 50);  //5000个线程块，维度是100*50
dim3 DimBlock(4, 8, 8);  //每个线层块内包含256个线程，线程块内的维度是4*8*8

怎么计算线程号？

• 使用N个线程块，每个线程块只有一个线程，即

dim3 dimGrid(N);
dim3 dimBlock(1);

此时的线程号计算方式就是

threadId = blockIdx.x;

其中，threadId的取值为0到N-1.对于这种情况，我们可以将其看作是一个列向量，列向量中的每一行对应一个线程块。列向量中每一行只有一个元素，对应一个线程。

• 使用MxN个线程块，每个线程块1个线程

由于线程块是2维的，故可以看作是一个M*N的2维矩阵，其线程号有两个维度，即：

dim3 dimGrid(M,N);
dim3 dimBlock(1);

其中

blockIdx.x 取值0到M-1
blockIdx.y 取值0到N-1

这种情况一般用于处理2维数据结构，比如2维图像。每个像素用一个线程来处理，此时需要线程号来映射图像像素的对应位置，如

pos = blockIdx.y * gridDim.x + blockIdx.x; //其中gridDim.x等于M

•  使用一个线程块，该线程块具有N个线程，即

dim3 dimGrid(1);
dim3 dimBlock(N);

此时线程号的计算方式维

threadId = threadIdx.x;

其中，threadId的范围是0到N-1，这种情况，可以看作是一个行向量，向量中的每一个元素的对应一个线程。

• 使用M个线程块，每个线程块内含有N个线程，即

dim3 dimGrid(M);
dim3 dimBlock(N);

这种情况，可以把它想象维二维矩阵，矩阵的行与线程块对应，矩阵的列于线程号对应，那么线程号的计算方式维为

threadId = threadIdx.x + blockIdx*blockDim.x;

上面其实就是把二维的索引空间转换为1为索引空间的过程。

•  使用MxN的二维线程块，每个线程块具有PxQ个线程，即

dim3 dimGrid(M, N);
dim3 dimBlock(P, Q);

这种情况是我们遇到最多的情况，特别适用于处理具有二维数据结构的算法，比如图像处理领域。

其索引有两个维度

threadId.x = blockIdx.x*blockDim.x+threadIdx.x;
threadId.y = blockIdx.y*blockDim.y+threadIdx.y;

上面的公式就是把线程和线程块的索引映射为图像像素坐标的计算方法。

6 CUDA应用例子

在掌握CUDA编程的基本语法，现在用一些小例子来真正入门CUDA编码。

6.1 查看我们GPU硬件配置

#include "device_launch_parameters.h"
#include <iostream>

int main()
{
    int deviceCount;
    cudaGetDeviceCount(&deviceCount);
    for(int i=0;i<deviceCount;i++)
    {
        cudaDeviceProp devProp;
        cudaGetDeviceProperties(&devProp, i);
        std::cout << "使用GPU device " << i << ": " << devProp.name << std::endl;
        std::cout << "设备全局内存总量： " << devProp.totalGlobalMem / 1024 / 1024 << "MB" << std::endl;
        std::cout << "SM的数量：" << devProp.multiProcessorCount << std::endl;
        std::cout << "每个线程块的共享内存大小：" << devProp.sharedMemPerBlock / 1024.0 << " KB" << std::endl;
        std::cout << "每个线程块的最大线程数：" << devProp.maxThreadsPerBlock << std::endl;
        std::cout << "设备上一个线程块（Block）种可用的32位寄存器数量： " << devProp.regsPerBlock << std::endl;
        std::cout << "每个EM的最大线程数：" << devProp.maxThreadsPerMultiProcessor << std::endl;
        std::cout << "每个EM的最大线程束数：" << devProp.maxThreadsPerMultiProcessor / 32 << std::endl;
        std::cout << "设备上多处理器的数量： " << devProp.multiProcessorCount << std::endl;
        std::cout << "======================================================" << std::endl;     
        
    }
    return 0;
}

 6.2 第一个计算任务

将两个元素数目为1024x1024的float数组相加。

用CPU串行编程思想完成这个任务：

 

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <stdlib.h>

#include <math.h>
#include <chrono>


int main()
{
    auto t1 = std::chrono::high_resolution_clock::now();
    float* A, * B, * C;
    int n = 1024 * 1024;
    int size = n * sizeof(float);
    A = (float*)malloc(size);
    B = (float*)malloc(size);
    C = (float*)malloc(size);

    for (int i = 0; i < n; i++)
    {
        A[i] = 90.0;
        B[i] = 10.0;
    }

    for (int i = 0; i < n; i++)
    {
        C[i] = A[i] + B[i];
    }

    float max_error = 0.0;
    for (int i = 0; i < n; i++)
    {
        max_error += fabs(100.0 - C[i]);
    }
    std::cout << "max_error is " << max_error << std::endl;
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << "total time is " << std::chrono::duration<double, std::milli>(t2 - t1).count() << " ms\n";
    return 0;
}

输出结果：

max_error is 0
total time is 28.9827 ms

 用GPU做并行计算：

 编码要点：

• 每个block中Thread数量最大不超过512
• 为了充分利用SM，Block数量尽可能多，>100

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <stdlib.h>

#include <math.h>
#include <chrono>

__global__ void Plus(float A[], float B[], float C[], int n)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    C[i] = A[i] + B[i];
}

int main()
{
    auto t1 = std::chrono::high_resolution_clock::now();
    float* A, * Ad, * B, * Bd, * C, * Cd;
    int n = 1024 * 1024;
    int size = n * sizeof(float);

    // CPU端分配内存
    A = (float*)malloc(size);
    B = (float*)malloc(size);
    C = (float*)malloc(size);

    // 初始化数组
    for (int i = 0; i < n; i++)
    {
        A[i] = 90.0;
        B[i] = 10.0;
    }

    // GPU端分配内存
    cudaMalloc((void**)&Ad, size);
    cudaMalloc((void**)&Bd, size);
    cudaMalloc((void**)&Cd, size);

    // CPU的数据拷贝到GPU端
    cudaMemcpy(Ad, A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);
    cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);

    // 定义kernel执行配置，（1024*1024/512）个block，每个block里面有512个线程
    dim3 dimBlock(512);
    dim3 dimGrid(n / 512);

    // 执行kernel
    Plus << <dimGrid, dimBlock >> > (Ad, Bd, Cd, n);

    // 将在GPU端计算好的结果拷贝回CPU端
    cudaMemcpy(C, Cd, size, cudaMemcpyDeviceToHost);

    // 校验误差
    float max_error = 0.0;
    for (int i = 0; i < n; i++)
    {
        max_error += fabs(100.0 - C[i]);
    }   

    // 释放CPU端、GPU端的内存
    free(A);
    free(B);
    free(C);
    cudaFree(Ad);
    cudaFree(Bd);
    cudaFree(Cd);
    std::cout << "max_error is " << max_error << std::endl;
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << "total time is " << std::chrono::duration<double, std::milli>(t2 - t1).count() << " ms\n";
    return 0;
}

输出：

max_error is 0
total time is 353.963 ms

从上面的例子可以看出，使用CUDA编程时，看不到for循环了。因为CPU编程的循环已经分散到各个thread上做了。从结果上看，CPU的循环计算速度比GPU计算块多了。原因在于CUDA中有大量的内存拷贝操作（数据传输花费了大量时间，而计算时间却非常少），如果计算量比较小的话，CPU计算会更合适一些。

6.3 第二个计算任务

矩阵的加法，即对两个矩阵对应的坐标元素相加后的结果存储到第三个矩阵对应的位置上。

CPU版本

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <stdlib.h>

#include <math.h>
#include <chrono>

#define ROWS 1024
#define COLS 1024

using namespace std;

int main()
{
    auto t1 = std::chrono::high_resolution_clock::now();
    int* A, ** A_ptr, * B, ** B_ptr, * C, ** C_ptr;
    int total_size = ROWS * COLS * sizeof(int);
    A = (int*)malloc(total_size);
    B = (int*)malloc(total_size);
    C = (int*)malloc(total_size);
    A_ptr = (int**)malloc(ROWS * sizeof(int*));
    B_ptr = (int**)malloc(ROWS * sizeof(int*));
    C_ptr = (int**)malloc(ROWS * sizeof(int*));

    //CPU一维数组初始化
    for (int i = 0; i < ROWS * COLS; i++)
    {
        A[i] = 80;
        B[i] = 20;
    }

    for (int i = 0; i < ROWS; i++)
    {
        A_ptr[i] = A + COLS * i;
        B_ptr[i] = B + COLS * i;
        C_ptr[i] = C + COLS * i;
    }

    for (int i = 0; i < ROWS; i++)
        for (int j = 0; j < COLS; j++)
        {
            C_ptr[i][j] = A_ptr[i][j] + B_ptr[i][j];
        }

    //检查结果
    int max_error = 0;
    for (int i = 0; i < ROWS * COLS; i++)
    {
        //cout << C[i] << endl;
        max_error += abs(100 - C[i]);
    }

    std::cout << "max_error is " << max_error << std::endl;
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << "total time is " << std::chrono::duration<double, std::milli>(t2 - t1).count() << " ms\n";
    return 0;
}

输出

max_error is 0
total time is 7.5455 ms

GPU版本

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <stdlib.h>

#include <math.h>
#include <chrono>

#define Row 1024
#define Col 1024

using namespace std;

__global__ void addKernel(int** C, int** A, int** B)
{
    int idx = threadIdx.x + blockDim.x * blockIdx.x;
    int idy = threadIdx.y + blockDim.y * blockIdx.y;
    if (idx < Col && idy < Row) {
        C[idy][idx] = A[idy][idx] + B[idy][idx];
    }
}

int main()
{
    auto t1 = std::chrono::high_resolution_clock::now();
    
    int** A = (int**)malloc(sizeof(int*) * Row);
    int** B = (int**)malloc(sizeof(int*) * Row);
    int** C = (int**)malloc(sizeof(int*) * Row);
    int* dataA = (int*)malloc(sizeof(int) * Row * Col);
    int* dataB = (int*)malloc(sizeof(int) * Row * Col);
    int* dataC = (int*)malloc(sizeof(int) * Row * Col);
    int** d_A;
    int** d_B;
    int** d_C;
    int* d_dataA;
    int* d_dataB;
    int* d_dataC;
    //malloc device memory
    cudaMalloc((void**)&d_A, sizeof(int**) * Row);
    cudaMalloc((void**)&d_B, sizeof(int**) * Row);
    cudaMalloc((void**)&d_C, sizeof(int**) * Row);
    cudaMalloc((void**)&d_dataA, sizeof(int) * Row * Col);
    cudaMalloc((void**)&d_dataB, sizeof(int) * Row * Col);
    cudaMalloc((void**)&d_dataC, sizeof(int) * Row * Col);
    //set value
    for (int i = 0; i < Row * Col; i++) {
        dataA[i] = 90;
        dataB[i] = 10;
    }
    //将主机指针A指向设备数据位置，目的是让设备二级指针能够指向设备数据一级指针
    //A 和  dataA 都传到了设备上，但是二者还没有建立对应关系
    for (int i = 0; i < Row; i++) {
        A[i] = d_dataA + Col * i;
        B[i] = d_dataB + Col * i;
        C[i] = d_dataC + Col * i;
    }

    cudaMemcpy(d_A, A, sizeof(int*) * Row, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, B, sizeof(int*) * Row, cudaMemcpyHostToDevice);
    cudaMemcpy(d_C, C, sizeof(int*) * Row, cudaMemcpyHostToDevice);
    cudaMemcpy(d_dataA, dataA, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    cudaMemcpy(d_dataB, dataB, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    dim3 threadPerBlock(16, 16);
    dim3 blockNumber((Col + threadPerBlock.x - 1) / threadPerBlock.x, (Row + threadPerBlock.y - 1) / threadPerBlock.y);
    printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);
    addKernel << <blockNumber, threadPerBlock >> > (d_C, d_A, d_B);
    //拷贝计算数据-一级数据指针
    cudaMemcpy(dataC, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);

    int max_error = 0;
    for (int i = 0; i < Row * Col; i++)
    {
        //printf("%d\n", dataC[i]);
        max_error += abs(100 - dataC[i]);
    }

    //释放内存
    free(A);
    free(B);
    free(C);
    free(dataA);
    free(dataB);
    free(dataC);
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    cudaFree(d_dataA);
    cudaFree(d_dataB);
    cudaFree(d_dataC);

    std::cout << "max_error is " << max_error << std::endl;
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << "total time is " << std::chrono::duration<double, std::milli>(t2 - t1).count() << " ms\n";
    return 0;
}

输出：

Block(16,16)   Grid(64,64).
max_error is 0
total time is 216.899 ms

 从结果上看，CPU计算时间还是比GPU的计算时间短。这里需要指出的是，这种二维数组的程序写法的效率并不高（虽然比较符合人的思维方式），因为我们做了两次内存访问操作。所以一般而言，做高性能计算一般不会采取这种编程方式。

6.4 第三个计算任务

计算一个更加复杂的任务，矩阵乘法。

回顾一下矩阵乘法：两个矩阵相乘，左矩阵第一行乘以右矩阵第一列（对应元素相乘），乘完之后相加，即为结果的第一行第一列的数，依次往下算，直到计算完成所有矩阵的元素。

 CPU版本

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <stdlib.h>

#include <math.h>
#include <chrono>

#define ROWS 1024
#define COLS 1024

using namespace std;

void matrix_mul_cpu(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.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()
{
    auto t1 = std::chrono::high_resolution_clock::now();
    
    float* A, * B, * C;
    int total_size = ROWS * COLS * sizeof(float);
    A = (float*)malloc(total_size);
    B = (float*)malloc(total_size);
    C = (float*)malloc(total_size);

    //CPU一维数组初始化
    for (int i = 0; i < ROWS * COLS; i++)
    {
        A[i] = 80.0;
        B[i] = 20.0;
    }

    matrix_mul_cpu(A, B, C, COLS);

   
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << "total time is " << std::chrono::duration<double, std::milli>(t2 - t1).count() << " ms\n";
    return 0;
}

输出：

total time is 3032.46 ms

梳理一下CUDA求解矩阵乘法的思路：因为C=AxB，我们利用每个线程求解C矩阵每个(x,y)的元素，每个线程载入A的一行和B的一列，变量各自行列元素吗，对A，B对应的元素做一次乘法和一次加法。

 GPU版本

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <stdlib.h>

#include <math.h>
#include <chrono>

#define Row 1024
#define Col 1024

using namespace std;

__global__ void matrix_mul_gpu(int* M, int* N, int* P, int width)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;

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

int main()
{
    auto t1 = std::chrono::high_resolution_clock::now();
    
    int* A = (int*)malloc(sizeof(int) * Row * Col);
    int* B = (int*)malloc(sizeof(int) * Row * Col);
    int* C = (int*)malloc(sizeof(int) * Row * Col);
    //malloc device memory
    int* d_dataA, * d_dataB, * d_dataC;
    cudaMalloc((void**)&d_dataA, sizeof(int) * Row * Col);
    cudaMalloc((void**)&d_dataB, sizeof(int) * Row * Col);
    cudaMalloc((void**)&d_dataC, sizeof(int) * Row * Col);
    //set value
    for (int i = 0; i < Row * Col; i++) {
        A[i] = 90;
        B[i] = 10;
    }

    cudaMemcpy(d_dataA, A, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    cudaMemcpy(d_dataB, B, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    dim3 threadPerBlock(16, 16);
    dim3 blockNumber((Col + threadPerBlock.x - 1) / threadPerBlock.x, (Row + threadPerBlock.y - 1) / threadPerBlock.y);
    printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);
    matrix_mul_gpu << <blockNumber, threadPerBlock >> > (d_dataA, d_dataB, d_dataC, Col);
    //拷贝计算数据-一级数据指针
    cudaMemcpy(C, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);

    //释放内存
    free(A);
    free(B);
    free(C);
    cudaFree(d_dataA);
    cudaFree(d_dataB);
    cudaFree(d_dataC);
   
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << "total time is " << std::chrono::duration<double, std::milli>(t2 - t1).count() << " ms\n";
    return 0;
}

输出：

Block(16,16)   Grid(64,64).
total time is 164.542 ms

从这个矩阵乘法任务可以看出，我们通过GPU进行并行计算的方式仅花了0.16秒，但CPU却花费了3秒，计算速度提升了近20倍。可见并行计算的威力。