CUDA编程——矩阵乘法的串行和两种并行实现

一.CUDA是什么

这里仅简单介绍一下主要概念,如下:
1.主机
将CPU及系统的内存(内存条)称为主机。
2.设备
将GPU及GPU本身的显示内存称为设备。
3.线程(Thread)
一般通过GPU的一个核进行处理。
4.线程块(Block)
1. 由多个线程组成。
2. 各block是并行执行的,block间无法通信,也没有执行顺序。
3. 注意线程块的数量限制为不超过65535(硬件限制)。
5.线程格(Grid)
由多个线程块组成。
1
6.线程束
在CUDA架构中,线程束是指一个包含32个线程的集合,这个线程集合被“编织在一起”并且“步调一致”的形式执行。在程序中的每一行,线程束中的每个线程都将在不同数据上执行相同的命令。
7.核函数(Kernel)
1. 在GPU上执行的函数通常称为核函数。
2. 一般通过标识符__global__修饰,调用通过<<<参数1,参数2>>>,用于说明内核函数中的线程数量,以及线程是如何组织的。
3. 以线程格(Grid)的形式组织,每个线程格由若干个线程块(block)组成,而每个线程块又由若干个线程(thread)组成。
4. 是以block为单位执行的。
5. 只能在主机端代码中调用。
6. 调用时必须声明内核函数的执行参数。
7. 在编程时,必须先为kernel函数中用到的数组或变量分配好足够的空间,再调用kernel函数,否则在GPU计算时会发生错误。
8.dim3结构类型
1. dim3是基于uint3定义的矢量类型,相当于由3个unsigned int型组成的结构体。uint3类型有三个数据成员unsigned int x; unsigned int y; unsigned int z;
2. 可使用于一维、二维或三维的索引来标识线程,构成一维、二维或三维线程块。
3. dim3结构类型变量用在核函数调用的<<<,>>>中。
4. 相关的几个内置变量
threadIdx,获取线程thread的ID索引;
blockIdx,线程块的ID索引;
blockDim,线程块的维度;
gridDim,线程格的维度。
5. 对于一维的block,线程的threadID=threadIdx.x。
6. 对于大小为(blockDim.x, blockDim.y)的 二维 block,线程的threadID=threadIdx.x+threadIdx.y*blockDim.x。
7. 对于大小为(blockDim.x, blockDim.y, blockDim.z)的 三维 block,线程的threadID=threadIdx.x+threadIdx.y*blockDim.x+threadIdx.z*blockDim.x*blockDim.y。
8. 对于计算线程索引偏移增量为已启动线程的总数。如stride = blockDim.x * gridDim.x; threadId += stride。
9.函数修饰符
1. __global__,表明被修饰的函数在设备上执行,但在主机上调用。
2. __device__,表明被修饰的函数在设备上执行,但只能在其他__device__函数或者__global__函数中调用。
10.同步方法__syncthreads()
确保线程块中的每个线程都执行完__syscthreads()前面的语句后,才会执行下一条语句。
注意:
1. 当线程块的数量为GPU中处理数量的2倍时,将达到最优性能。
2. 核函数执行的第一个计算就是计算输入数据的偏移。每个线程的起始偏移都是0到线程数量减1之间的某个值。然后,对偏移的增量为已启动线程的总数。

更具体的原理可参考这篇博文CUDA编程入门:向量加法和矩阵乘法

二.矩阵乘法算法设计思路

A的大小为[10*blocksize][10*blocksize], B的大小为[10*blocksize][20*blocksize],求C矩阵的结果。

1.串行矩阵乘法

void searial(int *A, int *B, int *C)
{
    for (int i = 0; i < 10 * BLOCK_SIZE; i++)
    {
        for (int j = 0; j < 20 * BLOCK_SIZE; j++)
        {
            int sum = 0;

            for (int k = 0; k < 10 * BLOCK_SIZE; k++)
            {
                sum += A[i * 10 * BLOCK_SIZE + k] * 
                             B[k * 20 * BLOCK_SIZE + j];
            }

            C[i * 20 * BLOCK_SIZE + j] = sum;
        }
    }
}

2.并行不分块计算矩阵乘法,具体方法为每个thread计算C矩阵中的一个元素,计算得到row值和col值后即可直接for循环计算

__global__
void deviceParallel1(int *A, int *B, int *C)
{

    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    int sum = 0;

    for (int i = 0; i < 10 * BLOCK_SIZE; i++)
    {
        sum += A[row * 10 * BLOCK_SIZE + i] * B[i * 20 * BLOCK_SIZE + col];
    }

    C[row * 20 * BLOCK_SIZE + col] = sum;
}

void parallel1(int *A, int *B, int *C)
{
    int *CA, *CB, *CC;

    cudaMalloc(&CA, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE);
    cudaMalloc(&CB, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE); 
    cudaMalloc(&CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);

    cudaMemcpy(CA, A, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE, 
                                        cudaMemcpyHostToDevice);

    cudaMemcpy(CB, B, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE, 
                                        cudaMemcpyHostToDevice);

    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE); 
    dim3 dimGrid(20, 10);

    deviceParallel1<<<dimBlock, dimGrid>>>(CA, CB, CC);

    cudaThreadSynchronize();

    cudaMemcpy(C, CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE, 
                                        cudaMemcpyDeviceToHost);

    cudaFree(CA);
    cudaFree(CB);
    cudaFree(CC);
}

3.并行分块计算矩阵乘法,具体方法为根据线程块号和块内线程号,循环遍历得到所有的子矩阵,在主要循环计算步骤前后添加同步线程语句,最后将结果写回到C矩阵

__global__
void deviceParallel2(int *A, int *B, int *C)
{
    //获得线程块号
    int blkRow = blockIdx.y; 
    int blkCol = blockIdx.x;

    //获得块内的线程号 
    int row = threadIdx.y; 
    int col = threadIdx.x;

    int var = 0;

    //循环,遍历所有子矩阵
    for (int i = 0; i < 10; i++) 
    {   
        const int *ASub = A + blkRow * BLOCK_SIZE * 10 + i * BLOCK_SIZE; 
        const int *BSub = B + i * BLOCK_SIZE * 20 + blkCol * BLOCK_SIZE;

        __shared__ int Ads[BLOCK_SIZE][BLOCK_SIZE]; 
        __shared__ int Bds[BLOCK_SIZE][BLOCK_SIZE];

        Ads[row][col] = *(ASub + row * BLOCK_SIZE * 10 + col); 
        Bds[row][col] = *(BSub + row * BLOCK_SIZE * 20 + col);

        __syncthreads();

        for (int i = 0; i < BLOCK_SIZE; i++) 
        {
            var += Ads[row][i] * Bds[i][col]; 
        }

        __syncthreads();
    }

    int *CSub = C + blkRow * BLOCK_SIZE * 20 + blkCol * BLOCK_SIZE;

    *(CSub + row * BLOCK_SIZE * 20 + col) = var;
}

void parallel2(int *A, int *B, int *C)
{
    int *CA, *CB, *CC;

    cudaMalloc(&CA, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE);
    cudaMalloc(&CB, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE); 
    cudaMalloc(&CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);

    cudaMemcpy(CA, A, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE, 
                                        cudaMemcpyHostToDevice);

    cudaMemcpy(CB, B, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE, 
                                        cudaMemcpyHostToDevice);

    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE); 
    dim3 dimGrid(20, 10);

    deviceParallel2<<<dimBlock, dimGrid>>>(CA, CB, CC);

    cudaThreadSynchronize();

    cudaMemcpy(C, CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE, 
                                        cudaMemcpyDeviceToHost);

    cudaFree(CA);
    cudaFree(CB);
    cudaFree(CC);
}

三.运行结果及加速比计算

2
计算加速比:(blocksize设置为32,若设置为16,则计算结果在时间上比较不明显,而64的话有可能会导致内存溢出)
第一种不分块并行算法:speedup = 656361/412700 = 1.59
第二种分块并行算法:speedup = 656361/5464 = 120.12

四.源代码

#include <stdio.h> 
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#include <iostream>
#include <cuda_runtime.h>

using namespace std;

#define BLOCK_SIZE 32

typedef void (*multiply)(int *A, int *B, int *C);

double getTime(int *A, int *B, int *C, multiply mul)
{
    timeval start, finish;

    gettimeofday(&start, 0);

    mul(A, B, C);

    gettimeofday(&finish, 0);

    double interval = 1e6 * (finish.tv_sec - start.tv_sec) + finish.tv_usec - start.tv_usec;

    return interval;
}

void searial(int *A, int *B, int *C)
{
    for (int i = 0; i < 10 * BLOCK_SIZE; i++)
    {
        for (int j = 0; j < 20 * BLOCK_SIZE; j++)
        {
            int sum = 0;

            for (int k = 0; k < 10 * BLOCK_SIZE; k++)
            {
                sum += A[i * 10 * BLOCK_SIZE + k] * 
                             B[k * 20 * BLOCK_SIZE + j];
            }

            C[i * 20 * BLOCK_SIZE + j] = sum;
        }
    }
}

__global__
void deviceParallel1(int *A, int *B, int *C)
{

    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    int sum = 0;

    for (int i = 0; i < 10 * BLOCK_SIZE; i++)
    {
        sum += A[row * 10 * BLOCK_SIZE + i] * B[i * 20 * BLOCK_SIZE + col];
    }

    C[row * 20 * BLOCK_SIZE + col] = sum;
}

void parallel1(int *A, int *B, int *C)
{
    int *CA, *CB, *CC;

    cudaMalloc(&CA, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE);
    cudaMalloc(&CB, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE); 
    cudaMalloc(&CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);

    cudaMemcpy(CA, A, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE, 
                                        cudaMemcpyHostToDevice);

    cudaMemcpy(CB, B, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE, 
                                        cudaMemcpyHostToDevice);

    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE); 
    dim3 dimGrid(20, 10);

    deviceParallel1<<<dimBlock, dimGrid>>>(CA, CB, CC);

    cudaThreadSynchronize();

    cudaMemcpy(C, CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE, 
                                        cudaMemcpyDeviceToHost);

    cudaFree(CA);
    cudaFree(CB);
    cudaFree(CC);
}

__global__
void deviceParallel2(int *A, int *B, int *C)
{
    //获得线程块号
    int blkRow = blockIdx.y; 
    int blkCol = blockIdx.x;

    //获得块内的线程号 
    int row = threadIdx.y; 
    int col = threadIdx.x;

    int var = 0;

    //循环,遍历所有子矩阵
    for (int i = 0; i < 10; i++) 
    {   
        const int *ASub = A + blkRow * BLOCK_SIZE * 10 + i * BLOCK_SIZE; 
        const int *BSub = B + i * BLOCK_SIZE * 20 + blkCol * BLOCK_SIZE;

        __shared__ int Ads[BLOCK_SIZE][BLOCK_SIZE]; 
        __shared__ int Bds[BLOCK_SIZE][BLOCK_SIZE];

        Ads[row][col] = *(ASub + row * BLOCK_SIZE * 10 + col); 
        Bds[row][col] = *(BSub + row * BLOCK_SIZE * 20 + col);

        __syncthreads();

        for (int i = 0; i < BLOCK_SIZE; i++) 
        {
            var += Ads[row][i] * Bds[i][col]; 
        }

        __syncthreads();
    }

    int *CSub = C + blkRow * BLOCK_SIZE * 20 + blkCol * BLOCK_SIZE;

    *(CSub + row * BLOCK_SIZE * 20 + col) = var;
}

void parallel2(int *A, int *B, int *C)
{
    int *CA, *CB, *CC;

    cudaMalloc(&CA, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE);
    cudaMalloc(&CB, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE); 
    cudaMalloc(&CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);

    cudaMemcpy(CA, A, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE, 
                                        cudaMemcpyHostToDevice);

    cudaMemcpy(CB, B, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE, 
                                        cudaMemcpyHostToDevice);

    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE); 
    dim3 dimGrid(20, 10);

    deviceParallel2<<<dimBlock, dimGrid>>>(CA, CB, CC);

    cudaThreadSynchronize();

    cudaMemcpy(C, CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE, 
                                        cudaMemcpyDeviceToHost);

    cudaFree(CA);
    cudaFree(CB);
    cudaFree(CC);
}

void read(int *M, int row, int col)
{
    srand((unsigned)time(NULL));

    for (int i = 0; i < row; i++)
    {
        for (int j = 0; j < col; j++)
        {
            M[i * col + j] = rand() % 1000;
        }
    }
}


int main(int argc, char const *argv[])
{

    int *A = new int[10 * BLOCK_SIZE * 10 * BLOCK_SIZE];

    if (A == NULL)
    {
        cerr << "can not malloc A" << endl; exit(1);
    }

    int *B = new int[10 * BLOCK_SIZE * 20 * BLOCK_SIZE];

    if (B == NULL)
    {
        cerr << "can not malloc B" << endl; exit(1);
    }

    int *C1 = new int[10 * BLOCK_SIZE * 20 * BLOCK_SIZE];

    if (C1 == NULL)
    {
        cerr << "can not malloc C" << endl; exit(1);
    }

    int *C2 = new int[10 * BLOCK_SIZE * 20 * BLOCK_SIZE];

    if (C2 == NULL)
    {
        cerr << "can not malloc C" << endl; exit(1);
    }

    int *C3 = new int[10 * BLOCK_SIZE * 20 * BLOCK_SIZE];

    if (C3 == NULL)
    {
        cerr << "can not malloc C" << endl; exit(1);
    }

    // 读取矩阵数据

    read(A, 10 * BLOCK_SIZE, 10 * BLOCK_SIZE);

    read(B, 10 * BLOCK_SIZE, 20 * BLOCK_SIZE);

    cout << "Serial Time = " << getTime(A, B, C1, searial) << " ps." << endl;

    cout << "Parallel1 Time = " << getTime(A, B, C2, parallel1) << " ps." << endl;

    cout << "Parallel2 Time = " << getTime(A, B, C3, parallel2) << " ps." << endl;

    delete[] A;
    delete[] B;
    delete[] C1;
    delete[] C2;
    delete[] C3;

    return 0;
}
已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页