GPU矩阵并行相乘

原文链接

实验介绍

1. 相对CPU来说，GPU更适合处理高度并行化的程序，此次实验借助CUDA架构，C++编码实现在GPU的矩阵快速相乘，
2. 实验中用到了CUDA的相关知识，如cudaMalloc，cudaMemcpy，cudaFree；clock_t，gettimeofday计算运行时间；
3. 线程块二维分布和一个线程块的线程数为256。
4. 与在CPU中的完成速度对比。
5. 采用内核函数，运用GPU的并行处理，对两个矩阵进行相乘（矩阵采用一维数组表示），矩阵采用随机函数rand()生成。

GPU简介

GPU结构

• NVIDIA的GPU在浮点运算能力上，吊打了Intel的CPU。其原因来自于CPU和GPU结构上的差异。
• 如下图所示，CPU仅仅具有有限的核心数量。
• 相比于GPU，CPU的核心属于“少而精”的存在，核心数虽然很少，
• 但是每个核心的性能很强，适合处理具有很多分支的复杂的逻辑。
• 近些年来，CPU中集成了一些并行指令集，如SSE、AVX等，其中AVX可以同时处理256位(32个字节)，
• 可以大大加速并行计算。但是相比于GPU，还是小巫见大巫。
  [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LptuRAXO-1626183587924)(https://i.ibb.co/4ssmDBs/1.jpg)]

CUDA简介

• CUDA(Compute Unified Device Architecture，计算统一设备架构)，竞争对手OpenCL(from 2008，苹果公司)。
• CUDA 是NVIDIA专有的，即只能用Nvidia的GPU。
• OpenCL是所有主流媒介采用的一直标准，OpenCL可以在所有平台(Nvidia, AMD等)执行，
• 但是否能具有好的运行效果会有差异，同一时刻CUDA更快，CUDA未来会比OpenCL发展更快。

线程讲解

• CUDA编程是一个多线程编程，数个线程(Thread)组成一个线程块(Block)，所有线程块组成一个线程网格(Grid)，
• 图中的线程块，以及线程块中的线程，是按照2维的方式排布的。
• 实际上，CUDA编程模型允许使用1维、2维、3维三种方式来排布。
• 另外，即使线程块使用的是1维排布，线程块中的线程也不一定要按照1维排，而是可以任意排布。
  [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ndY1Ewv4-1626183587949)(https://i.ibb.co/qWvGQDT/2.jpg)]
• 目前的GPU限制一个线程块中，最多可以安排1024个线程。
• 一个线程块用多少线程，以及一个线程网格用多少线程块，是程序员可以自由安排的。
• 一般线程块中线程的数量被安排为32的倍数，选用256是比较合适的。
• 在线程数定下来之后，一般根据数据的排布情况来确定线程块的个数。（1维排列256，2维排列（16，16））
• 例如：一个数组的长度为4096，安排每个线程处理一个元素。
• 如果安排一个线程块为256个线程，则需要4096/256=16个线程块。

内核函数

内核函数是CUDA 每个线程 执行的函数，它运行在GPU设备上。CUDA使用扩展的C语言编写内核函数，关键字为__global__。内核函数返回值只能是void。

• 定义格式：__global__ void 函数名（参数……）{ 程序指令集合 }
• 主函数调用的格式：函数名<<<blocksPerGrid, threadsPerBock>>>（参数……）
• blocksPerGrid：每个网格中进程块的排布方式（可以采用1维或2维）
• threadsPerBock:每个进程块中进程的排布方式（可以采用1维或2维）
  内核函数举例

_global void VecAdd(double a[][],double b[][],double c[][]){
    int x = blockIdx.x * blockDim.x  + threadIdx.x;//当前列址
    int y = blockIdy.y * blockDim.y + threadIdy.y;//当前行址0099
    if(i < N&& y < N){
        c[j][i] = a[j][i] + b[j][i]
    }
}
int main(){
    dim3 threadsPerBlock(16,16);//每个线程块内部排布
    dim3 blocksPerGrid(N / threadsPerBlock.x ,N / threadsPerBlock.y);//线程排布
    VecAdd<<<blocksPerGrid,threadsPerBlock>>>(A,B,C);
}

编程接口

• 使用NVCC编译CUDA程序
• CUDA程序使用NVCC编译器。NVCC提供了简单方便的接口，能够很好的同时处理主机端和设备端代码。
• 编译程序的命令：nvcc filename.cu –o filename

cuda主要函数

1. cudaMalloc

   cudaMalloc (void **devPtr, size_t size )

2. cudaMemcpy

   • 主机到设备：cudaMemcpy(d_A,h_A,nBytes,cudaMemcpyHostToDevice)

   • 设备到主机：cudaMemcpy(h_A,d_A,nBytes,cudaMemcpyDeviceToHost)

实践作业

• 编写一个矩阵乘法的GPU并行程序，
• 并且与对应规模的串行程序进行运行时间的比对（n=500，1000，1500，2000，3000，5000），
• 画出规模和时间对比图。
• 矩阵A（n，n） 矩阵B（n，n） C = A x B

内核函数

#include <iostream>

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include "cuda_runtime.h"
#include <cuda.h>
#include <sys/time.h> 

#include "device_launch_parameters.h"

#define thread_num 256//一个线程块的线程数
using namespace std;


const int N = 6000;//数组维数

const int blocks_num = (N + thread_num - 1) / thread_num;//线程块数

__global__ void mextix(int *da,int *db,int *dc)
{
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    int col = blockIdx.y * blockDim.y + threadIdx.y;
    if(row < N && col < N){
        dc[row*N+col] = 0;
        for(int i = 0;i < N;i++){
            dc[row*N+col] += da[row*N+i] * db[i*N+col];
        }
        
    }
   
}

//随机生成矩阵
void rands(int *a)
{
   for(int i = 0;i < N;i++){
       for(int j = 0;j < N;j++){
           a[i*N+j] = rand() % 10 + 1 ;
       }
   }
}
int main()
{
    int *a,*b,*c;
    int *da,*db,*dc;
    int size = N*N*sizeof(int);
    //freopen("out.txt","w",stdout);
    //分配空间
    a = (int*)malloc(size);
    b = (int*)malloc(size);
    c = (int*)malloc(size);

    //生成随机数组
    rands(a);
    rands(b);

    //分配内存 GPU申请空间所需时间
    clock_t t1 = clock();
    cudaMalloc((void**)&da,size);
    cudaMalloc((void**)&db,size);
    cudaMalloc((void**)&dc,size);
    //cudaMalloc((void**)&time,blocks_num*sizeof(clock_t)*2);
    clock_t t2 = clock();
    double ts = (double)(t2-t1);

    //CLOCKS_PER_SEC表示一秒钟内CPU运行的时钟周期数
    printf("GPU divide costtime : %lf ms\n",ts/CLOCKS_PER_SEC*1000);
    
    //存到GPU
    cudaMemcpy(da,a,size,cudaMemcpyHostToDevice);
    cudaMemcpy(db,b,size,cudaMemcpyHostToDevice);

    /*
        GPU运算  并行运算时间
        计算代码运行时间
    */
    timeval start,finish1,finish2;
    gettimeofday(&start,0);//获得当前精确时间
    dim3 dg(16,16);
    dim3 dbs((N+dg.x-1)/dg.x,(N+dg.y-1)/dg.y);
    gettimeofday(&finish1,0);
    mextix<<<dbs,dg>>>(da,db,dc);
    gettimeofday(&finish2, 0);//获得当前精确时间
    double cost1 = 1e6 * (finish2.tv_sec - start.tv_sec) + finish2.tv_usec - start.tv_usec;//微秒
    double cs = 1e6*(finish1.tv_sec - start.tv_sec) + (finish1.tv_usec - start.tv_usec);
    /*  
        timeval
        {
            time_t tv_sec;  //秒 [long int]
            suseconds_t tv_usec;  //微秒 [long int]
        };
   */

    //从GPU取回
    cudaMemcpy(c,dc,size,cudaMemcpyDeviceToHost);

    //GPU运算时间
    printf("GPUCost time : %lf ms\n",cost1/1e3);
    printf("GPU divdided time : %lf ms\n",cs/1e3);
    // printf("GPUAnswer : \n");
    // for(int i = 0;i < N;i++){
    //     for(int j = 0;j < N;j++){
    //         printf("%d ",c[i*N+j]);
    //         //printf("1");
    //     }
    //     printf("\n");
    // }
    //释放内存
    cudaFree(da);
    cudaFree(db);
    cudaFree(dc);
    CPU计算
    clock_t st = clock();
    for(int i = 0;i < N;i++){
        for(int j = 0;j < N;j++){
            c[i*N+j] = 0;
            for(int k = 0;k < N;k++){
                c[i*N+j] += a[i*N+k] * b[k*N+j];
            }
            
        }
        
    }
    clock_t ed = clock();
    double ends = (double)(ed-st);
    CPU运算时间
    printf("CPUCost time : %lf ms\n",ends/CLOCKS_PER_SEC*1000);
    // printf("CPUAnswer : \n");
    // for(int i = 0;i < N;i++){
    //     for(int j = 0;j < N;j++){
    //         printf("%d ",c[i*N+j]);
    //     }
    //     printf("\n");
    // }
    return 0;
}

性能对比分析

• GPU运行时间与n呈线性关系，运行时间随n的增大而增大；
• CPU运行时间与n呈指数关系，运行时间随n的增大而增大。
• GPU运行时间在毫秒级，而CPU则在秒级，GPU运行时间远远小于CPU。

分析可得，因为GPU采用线程并行处理矩阵相乘，而CPU采用串行一个个依次算，所以GPU运行时间会更短。