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


博主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分散归约向量内积

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

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

#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)
  • 5
    点赞
  • 24
    收藏
    觉得还不错? 一键收藏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值