CUDA学习(十)：向量内积的多种方法实现

最新推荐文章于 2024-09-26 16:53:57 发布

hjxu2016

最新推荐文章于 2024-09-26 16:53:57 发布

阅读量2.5k

点赞数 5

分类专栏：编程语言|CUDA入门

本文链接：https://blog.csdn.net/hjxu2016/article/details/109781166

版权

编程语言|CUDA入门专栏收录该内容

19 篇文章 149 订阅

订阅专栏

博主CUDA学习系列汇总传送门（持续更新）：编程语言|CUDA入门

文章目录

本文章为《 GPU编程与优化大众高性能计算》的读书笔记，例子也都取自书中教程。

向量内积运算中各元素的运算间存在简单联系，即需要累加所有元素乘积结果。
向量内积结果的累加过程称为归约（reduction）
向量内积运算在数学上的计算公式为
在这里插入图片描述也就是A和B的向量长度相等，相同索引的元素相乘，然后将乘积结果累加，得到结果C

一、CPU上实现向量内积

#include <stdio.h>
#include <iostream>

const int N = 2048;
/* 一、cpu 向量内积 */
template <typename T>
void dot_cpu(T *a, T *b, T *c, int n)
{
    double dTemp = 0;
    for (int i=0; i<n; ++i)
    {
        dTemp += a[i] * b[i];
    }
    *c = dTemp;
}

int main()
{
    float a[N], b[N];
    float c = 0;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i * 1.0;
        b[i] = 1.0;
    }
    dot_cpu(a, b, &c, N);
    printf("a * b = %f\n", c);
    printf("hello world\n");
}

二、GPU下单Block分散归约向量内积

各线程独立完成所有向量元素的乘积运算，然后将乘积的结果累加到向量内积结果。
thread是私有的，block内的threads要进行数据通信，必须使用共享内存或者全局内存。
共享内存的访存速度比全局快一个量级，所以我们这边采用共享内存。

假设单Block 8 线程求内积， 0-7分别表示各线程内内积结果
在这里插入图片描述该方法存在缺陷

单线程访问全局的相邻两个单元，违背访问对齐原则，关于对齐可参考gpu显存（全局内存）在使用时数据对齐的问题
共享存储的访问更容易引起 back conflict
什么是back conflict？关于共享存储器，如果同一个warp的不同线程访问到同一个bank中的数据，就会发生 bank conflict。更多参考CUDA-Memory(存储)和bank-conflict

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

const int N = 2048;

const int threadnum = 32;//开32个线程

/* 二、单block 分散归约 */
template <typename T>
__global__ void dot_gpu_1(T *a, T *b, T *c, int n)
{
    __shared__ T tmp[threadnum];
    const int nThreadIdX = threadIdx.x; //线程ID索引号
    const int nBlockDimX = blockDim.x; // 一个block内开启的线程总数
    int nTid = nThreadIdX;
    double dTemp = 0.0;
    while (nTid < n)
    {
        dTemp += a[nTid] * b[nTid];
        nTid += nBlockDimX;
    }
    tmp[nThreadIdX] = dTemp; // 将每个线程中的内积放入到共享内存中
    __syncthreads(); // 同步操作，即等所有线程内上面的操作都执行完

    int i=2, j=1;
    while(i <= threadnum)
    {
        if (nThreadIdX%i == 0)
        {
            tmp[nThreadIdX] += tmp[nThreadIdX + j];
        }
        __syncthreads();
        i *= 2;
        j *= 2;
    }
    if (0 == nThreadIdX)
    {
        c[0] = tmp[0];
    }
}


int main()
{
    float a[N], b[N];
    float c = 0;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i * 1.0;
        b[i] = 1.0;
    }
    float *d_a = NULL, *d_b = NULL, *d_c = NULL;
    cudaMalloc(&d_a, N *sizeof(float));
    cudaMemcpy(d_a, a, N *sizeof(float), cudaMemcpyHostToDevice);

    cudaMalloc(&d_b, N *sizeof(float));
    cudaMemcpy(d_b, b, N *sizeof(float), cudaMemcpyHostToDevice);

    cudaMalloc(&d_c, sizeof(float));
    dot_gpu_1<<<1, threadnum>>>(d_a, d_b, d_c, N);
    cudaMemcpy(&c, d_c, sizeof(float), cudaMemcpyDeviceToHost);
        std::cout << c << std::endl;
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    printf("hello world\n");
    return 0;
}

三、单Block低线程归约向量内积

在这里插入图片描述在归约时，访问共享存储区时，不同线程尽可能少的访问同一个bank中的数据地址。
一块共享存储取中存在多个bank，访问共享内存中的数据时，尽可能离远点。

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

const int N = 2048;

const int threadnum = 32;//开32个线程


/* 三、单block 低线程归约向量内积*/
template <typename T>
__global__ void dot_gpu_2(T *a, T *b, T *c, int n)
{
    __shared__ T tmp[threadnum];
    const int nThreadIdX = threadIdx.x; //线程ID索引号
    const int nBlockDimX = blockDim.x; // 一个block内开启的线程总数
    int nTid = nThreadIdX;
    double dTemp = 0.0;
    while (nTid < n)
    {
        dTemp += a[nTid] * b[nTid];
        nTid += nBlockDimX;
    }
    tmp[nThreadIdX] = dTemp; // 将每个线程中的内积放入到共享内存中
    __syncthreads(); // 同步操作，即等所有线程内上面的操作都执行完

    int i = threadnum / 2;
    while(i != 0)
    {
        if (nThreadIdX < i)
        {
            tmp[nThreadIdX] += tmp[nThreadIdX + i];
        }
        __syncthreads();// 同步操作，即等所有线程内上面的操作都执行完
        i /= 2;
    }
    if (0 == nThreadIdX)
    {
        c[0] = tmp[0];
    }
}

int main()
{
    float a[N], b[N];
    float c = 0;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i * 1.0;
        b[i] = 1.0;
    }
    float *d_a = NULL, *d_b = NULL, *d_c = NULL;
    cudaMalloc(&d_a, N *sizeof(float));
    cudaMemcpy(d_a, a, N *sizeof(float), cudaMemcpyHostToDevice);

    cudaMalloc(&d_b, N *sizeof(float));
    cudaMemcpy(d_b, b, N *sizeof(float), cudaMemcpyHostToDevice);

    cudaMalloc(&d_c, sizeof(float));
    dot_gpu_2<<<1, threadnum>>>(d_a, d_b, d_c, N);
    cudaMemcpy(&c, d_c, sizeof(float), cudaMemcpyDeviceToHost);
        std::cout << c << std::endl;
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    printf("hello world\n");
    return 0;
}

四、多block向量内积

只有一个block，显然无法发挥GPU的真实性能。
这一节使用多block，block内部采用低线程归约的方法，每个block计算得到一个中间结果并组成一个数组，然后再将中间结果进行二次归约，二次归约可以在cpu上实现，也可以在Gpu上实现，这里只记录在cpu上实现的方法。

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

const int N = 2048;

const int nBlockNum = 16;//开32个block
const int threadnum = 32;//开32个线程


/* 三、单block 低线程归约向量内积*/

/* 四、多block多线程向量内积 */
template <typename T>
__global__ void dot_gpu_3(T *a, T *b, T *c, int n)
{
    __shared__ T aTmp[threadnum];
    const int nThreadIdX = threadIdx.x; //线程ID索引号
    const int nStep = gridDim.x * blockDim.x; // 跳步的步长，即所有线程的数量
    int nTidIdx = blockIdx.x * blockDim.x + threadIdx.x; // 当前线程在全局线程的索引

    double dTemp = 0.0;
    while (nTidIdx < n)
    {
        dTemp += a[nTidIdx] * b[nTidIdx];
        nTidIdx += nStep;
    }
    aTmp[nThreadIdX] = dTemp; // 将每个线程中的内积放入到对应block的共享内存中
    __syncthreads(); // 同步操作，即等所有线程内上面的操作都执行完

    int i = threadnum / 2;
    while (i != 0)
    {
        if(nThreadIdX < i)
        {
            aTmp[nThreadIdX] += aTmp[nThreadIdX + i];
        }
        __syncthreads(); // 同步操作，即等所有线程内上面的操作都执行完
        i /= 2;
    }

    if (0 == nThreadIdX)
    {
        c[blockIdx.x] = aTmp[0];
    }

}

int main()
{
     float a[N], b[N];
    float c = 0;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i * 1.0;
        b[i] = 1.0;
    }
    float *d_a = NULL, *d_b = NULL, *d_c = NULL;
    cudaMalloc(&d_a, N *sizeof(float));
    cudaMemcpy(d_a, a, N *sizeof(float), cudaMemcpyHostToDevice);

    cudaMalloc(&d_b, N *sizeof(float));
    cudaMemcpy(d_b, b, N *sizeof(float), cudaMemcpyHostToDevice);


    cudaMalloc(&d_c, sizeof(float) * nBlockNum);
    dot_gpu_3<<< nBlockNum, threadnum >>>(d_a, d_b, d_c, N);
    float c_tmp[nBlockNum];
    cudaMemcpy(&c_tmp, d_c, nBlockNum * sizeof(float), cudaMemcpyDeviceToHost);
    for (int i=0; i < nBlockNum; ++i)
    {
        c += c_tmp[i];
    }
 	cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    std::cout << c << std::endl;
    printf("hello world\n");
    return 0;
}

五、cublas库实现向量内积

/* 五、cublas 实现向量内积 */
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <iostream>
#include "cublas_v2.h"
#include "common.h"
const int N = 2048;


int main()
{
    float a[N], b[N];
    float c = 0.0;
    for(int i = 0; i < N; ++i) // 为数组a、b赋值
    {
        a[i] = i * 1.0;
        b[i] = 1.0;
    }
    float *d_a = NULL, *d_b = NULL;

    cublasHandle_t handle;
    cublasCreate_v2(&handle); // 创建句柄
    cudaMalloc(&d_a, N *sizeof(float));
    cublasSetVector(N, sizeof(float), a, 1, d_a, 1);

    cudaMalloc(&d_b, N *sizeof(float));
    cublasSetVector(N, sizeof(float), b, 1, d_b, 1);
    
    CHECK(cublasSdot_v2(handle, N, d_a, 1, d_b, 1, &c));

  
    CHECK(cudaFree(d_a));
    CHECK(cudaFree(d_b));
    cublasDestroy(handle); // 销毁句柄

    printf("a * b = %f\n", c);
    printf("hello world\n");
    return 0;
}

CHECK函数的定义如下

#define CHECK(status)                                                                                                  \
    do                                                                                                                 \
    {                                                                                                                  \
        auto ret = (status);                                                                                           \
        if (ret != 0)                                                                                                  \
        {                                                                                                              \
            std::cout << "Cuda failure: " << ret << std::endl;                                                 \
            abort();                                                                                                   \
        }                                                                                                              \
    } while (0)