【高性能计算】CUDA

可乐大牛

已于 2022-10-24 21:02:42 修改

阅读量700

点赞数 1

分类专栏：杂谈文章标签：其他

于 2022-09-08 10:57:12 首次发布

本文链接：https://blog.csdn.net/qq_44173974/article/details/126757788

版权

杂谈专栏收录该内容

38 篇文章 10 订阅

订阅专栏

学习来源

大佬链接1

原理

GPU中的几个基本概念
 GPU硬件结构
 thread/block的概念(grid和block不同的划分对应的ThreadID计算)
共享内存与全局内存

原理+实践

编程语言|CUDA入门

案例

注：

c++一般是cpp文件，用gcc编译；而cuda的程序一般是cu文件，用nvcc编译。其实cuda程序可以是cpp文件，但是处理起来会麻烦，这里就直接用cu文件了。

案例一

CUDA中的hello world！
helllo.cu，运行：nvcc hello.cu -o hello && ./hello

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


__global__ void test(void)  // __global__ 为CUDA的关键字,表示代码在设备端(GPU端)运行, 可以在CPU端被调用
{
	printf("Hello CUDA!\n");
}

int main()
{
	test <<<1, 1 >>> ();  // 函数调用,  <<< >>>中的第一个参数表示块的个数, 第二个参数表示每个线程块中线程的个数
						  // 这里是使用一个线程块,这个线程块中只有一个线程执行这个函数.
	cudaDeviceSynchronize(); // 会阻塞当前程序的执行，直到所有任务都处理完毕（这里的任务其实就是指的是所有的线程都已经执行完了kernel function）。
							 // 通俗讲,就是等待设备端的线程执行完成
							 // 一个线程块中可以有多个线程,GPU的线程是GPU的最小操作单元
							 
 	return 0;
}

案例二

CUAD程序设计流程的演示：

获取设备
分配显存
数据传输（从CPU到GPU）
核函数
数据传输（从GPU到CPU）
释放显存空间
重置设备（可以省略）

add.cu，运行：nvcc add.cu -o add && ./add

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


__global__ void add(int a, int b, int *c)
{
    *c = a + b;
}
int main() {
    int *c;
    int *dev_c;
    cudaError_t cudaStatus;
    cudaStatus = cudaMalloc(&dev_c, sizeof(int)); // 分配显存
    if(cudaSuccess != cudaStatus)
    {
        fprintf(stderr, "cuda melloc error!");
        return -1;
    }
    add<<<1, 1>>>(2, 7, dev_c);  // 核函数计算
    cudaStatus = cudaMemcpy(&c, dev_c, sizeof(int), cudaMemcpyDeviceToHost); //  数据传输（从GPU到CPU）
    printf("2 + 7 = %d \n", c);
    cudaFree(dev_c); // 释放显存空间
    return 0;
}

案例三

grid划分成1维，block划分为1维的情况下实现向量加法（渐进式的实现）
vector_add.cu，运行：nvcc vector_add.cu -o vector_add && ./vector_add

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

#define  N 10


/**======================================================================
1、单block单thread向量加 
一个kernel对应一个grid，一个grid包含一个block，一个block包含一个thread
=======================================================================**/
__global__ void vector_add_gpu_1(int *a, int *b, int *c, int n)
{
    for(int i=0; i<n; ++i)
    {
        c[i] = a[i] + b[i];
    }
}

/**====================================================
2、单block多thread向量加
每个block中的thread是并行的，当前设置的是4个，线程索引号分别为0,1,2,3
所以最开始计算c[0],c[1],c[2],c[3],然后每个线程跳步计算c[4],c[5],c[6],c[7], 以此类推
=====================================================**/
__global__ void vector_add_gpu_2(int *a, int *b, int *c, int n)
{
    int tid = threadIdx.x; //线程索引号
    printf("%d\t",tid);
    const int t_n = blockDim.x; // 一个block内的线程总数
    while (tid < n)
    {
        c[tid] = a[tid] + b[tid];
        tid += t_n;
    }
}


/**====================================================
3、多block多thread向量加 
其实和第二个类似，这个是标准的 grid划分成1维，block划分为1维的情况
=====================================================**/
__global__ void vector_add_gpu_3(int *a, int *b, int *c, int n)
{
    int tid = blockIdx.x * blockDim.x + threadIdx.x; // 获取线程索引
    const int t_n = gridDim.x * blockDim.x; // 跳步的步长，所有线程的数量
    while (tid < n)
    {
        c[tid] = a[tid] + b[tid];
        tid += t_n;
    }
}




int main() {
    int a[N], b[N], c[N];
    int *dev_a, *dev_b, *dev_c;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i;
        b[i] = i * i;
    }
    cudaMalloc(&dev_a, sizeof(int) * N);
    cudaMemcpy(dev_a, a, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_b, sizeof(int) * N);
    cudaMemcpy(dev_b, b, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_c, sizeof(int) * N);
    cudaMemcpy(dev_c, c, sizeof(int) * N, cudaMemcpyHostToDevice);

    //  <<<n,m>>> n表示grid中有多少个
//   vector_add_gpu_1<<<1, 1>>>(dev_a, dev_b, dev_c, N);
//    vector_add_gpu_2<<<1, 4>>>(dev_a, dev_b, dev_c, N);
    vector_add_gpu_3<<<2, 4>>>(dev_a, dev_b, dev_c, N);
    cudaMemcpy(c, dev_c, sizeof(int) * N, cudaMemcpyDeviceToHost);

    for(int i=0; i<N; ++i)
    {
        printf("%d + %d = %d \n", a[i], b[i], c[i]);
    }

    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);
    return 0;
}

案例四

从之前的案例三进一步出发，得到向量点乘操作的实现，就是向量逐元素相乘再相加（渐进式的实现）
vector_dot.cu，运行：nvcc vector_dot.cu -o vector_dot&& ./vector_dot

我们之前的向量加法是每个线程负责若干个元素相加的工作，每一次相加完成之后就往 global memory写结果，彼此之间互相独立；
我们本次的向量点积需要每个线程负责若干个元素相乘的工作，每一次相乘完成之后可以写出，但是因为最终需要将这些结果一起相加的，所以我们不写出，而是将这些结果累加保存，也就是说每个线程维护了一个临时的内积，那么我们最后将这些临时内积加起来就是最终的向量点积结果了，也就是需要进行线程通信，而要进行线程通信就必须使用共享内存或者全局内存了，考虑到速度，我们使用共享内存。

而临时结果的合并方式参考下图：
在这里插入图片描述

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


#define  N 8
#define BLOCK_NUM 1
#define THREAD_NUM 8



__global__ void vector_dot_gpu(int *a, int *b, int *c, int n)
{
    __shared__  int tmp[THREAD_NUM]; //每个block的共享内存
    const int tid = threadIdx.x; //当前block内的线程索引号
    const int t_n = blockDim.x; // 一个grid内的线程总数 当前只有一个block

    int global_tid=tid; //当前线程在整个grid中的索引号
    int t_sum=0; //当前线程完成的临时内积

    while (global_tid < n)
    {
        t_sum = a[global_tid] * b[global_tid];
        global_tid += t_n;
    }
    tmp[tid]=t_sum; //当前线程完成的临时内积 存到共享内存中
    __syncthreads(); // 同步操作，即等所有线程内上面的操作都执行完

    //分散归约过程
    int i=2, j=1;
    while(i <= THREAD_NUM)
    {
        if (tid%i == 0)
        {
            tmp[tid] += tmp[tid + j];
        }
        __syncthreads();
        i *= 2;
        j *= 2;
    }
    if (0 == tid)
    {
        c[0] = tmp[0];
    }
}



int main() {
    int a[N], b[N], c;
    int *dev_a, *dev_b, *dev_c;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i;
        b[i] = i;
    }
    cudaMalloc(&dev_a, sizeof(int) * N);
    cudaMemcpy(dev_a, a, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_b, sizeof(int) * N);
    cudaMemcpy(dev_b, b, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_c, sizeof(int));

    vector_dot_gpu<<<BLOCK_NUM, THREAD_NUM>>>(dev_a, dev_b, dev_c, N);
    cudaMemcpy(&c, dev_c, sizeof(int), cudaMemcpyDeviceToHost);

    printf("c=%d\n",c);

    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);
    return 0;
}

临时结果的合并方式（改）参考下图：
在这里插入图片描述

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

多block版本，每个block有一个临时内积，我们再度规约这些临时内积即可，这个规约的过程可以在GPU上实现，也可以在cpu上实现，下面是cpu上实现。

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


#define  N 8
#define BLOCK_NUM 2
#define THREAD_NUM 4



__global__ void vector_dot_gpu(int *a, int *b, int *c, int n)
{
    __shared__  int tmp[THREAD_NUM]; //每个block的共享内存
    const int tid = threadIdx.x; //当前block内的线程索引号
    const int t_n = gridDim.x*blockDim.x; // 一个grid内的线程总数

    int global_tid=blockIdx.x*blockDim.x+tid; //当前线程在整个grid中的索引号
    int t_sum=0; //当前线程完成的临时内积

    while (global_tid < n)
    {
        t_sum = a[global_tid] * b[global_tid];
        global_tid += t_n;
    }
    tmp[tid]=t_sum; //当前线程完成的临时内积 存到共享内存中
    __syncthreads(); // 同步操作，即等所有线程内上面的操作都执行完

    //低线程归约过程
    int i = THREAD_NUM / 2;
    while(i != 0)
    {
        if (tid < i)
        {
            tmp[tid] += tmp[tid + i];
        }
        __syncthreads();// 同步操作，即等所有线程内上面的操作都执行完
        i /= 2;
    }
    if (0 == tid)
    {
        c[blockIdx.x] = tmp[0];
    }
}



int main() {
    int a[N], b[N], c=0;
    int *dev_a, *dev_b, *dev_c;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i;
        b[i] = i;
    }
    cudaMalloc(&dev_a, sizeof(int) * N);
    cudaMemcpy(dev_a, a, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_b, sizeof(int) * N);
    cudaMemcpy(dev_b, b, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_c, sizeof(int)*BLOCK_NUM);

    vector_dot_gpu<<<BLOCK_NUM, THREAD_NUM>>>(dev_a, dev_b, dev_c, N);

    int c_tmp[BLOCK_NUM]={0};
    cudaMemcpy(&c_tmp, dev_c, sizeof(int)*BLOCK_NUM, cudaMemcpyDeviceToHost);
    for(int i=0;i<BLOCK_NUM;i++)
    {
        printf("%d\n",c_tmp[i]);
        c+=c_tmp[i];
    }    
    printf("c=%d\n",c);

    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);
    return 0;
}

GPU上实现二次规约过程


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


#define  N 8
#define BLOCK_NUM 2
#define THREAD_NUM 4



__global__ void vector_dot_gpu(int *a, int *b, int *c, int n)
{
    __shared__  int tmp[THREAD_NUM]; //每个block的共享内存
    const int tid = threadIdx.x; //当前block内的线程索引号
    const int t_n = gridDim.x*blockDim.x; // 一个grid内的线程总数

    int global_tid=blockIdx.x*blockDim.x+tid; //当前线程在整个grid中的索引号
    int t_sum=0; //当前线程完成的临时内积

    while (global_tid < n)
    {
        t_sum = a[global_tid] * b[global_tid];
        global_tid += t_n;
    }
    tmp[tid]=t_sum; //当前线程完成的临时内积 存到共享内存中
    __syncthreads(); // 同步操作，即等所有线程内上面的操作都执行完

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



int main() {
    int a[N], b[N], c;
    int *dev_a, *dev_b, *dev_c;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i;
        b[i] = i;
    }
    cudaMalloc(&dev_a, sizeof(int) * N);
    cudaMemcpy(dev_a, a, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_b, sizeof(int) * N);
    cudaMemcpy(dev_b, b, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_c, sizeof(int));
    cudaMemset(&dev_c, 0, sizeof(int));

    vector_dot_gpu<<<BLOCK_NUM, THREAD_NUM>>>(dev_a, dev_b, dev_c, N);

    cudaMemcpy(&c, dev_c, sizeof(int), cudaMemcpyDeviceToHost);  
    printf("c=%d\n",c);

    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);
    return 0;
}

注意到，我们上面block之间临时内积的累加采用的是c[0] += tmp[0]，其实这是存在风险的，多个线程并发访问c[0]，所以我们可以将它改成CUDA的原子函数实现atomicAdd(c, tmp[0]);，原子函数顾名思义就是线程安全的原子计算，那么话也说回来了，如果真的要用这个原子函数，我们其实没必要像上面那么麻烦了，直接就将每个线程的临时内积使用原子函数累加到c[0]中就好了，不需要共享内存，不需要规约了。注意：写法是简单了，但是原子操作对数据访问是串行的，频繁的原子操作会非常耗时。
也就是说，这边提到了两种思路，第一种是直接Block内归约，block间原子操作，第二种是直接使用原子操作代替两个规约操作，后者简单，但是速度慢，前者稍微简单一点，速度也均衡点。

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


#define  N 8
#define BLOCK_NUM 2
#define THREAD_NUM 4



__global__ void vector_dot_gpu(int *a, int *b, int *c, int n)
{
    const int tid = threadIdx.x; //当前block内的线程索引号
    const int t_n = gridDim.x*blockDim.x; // 一个grid内的线程总数

    int global_tid=blockIdx.x*blockDim.x+tid; //当前线程在整个grid中的索引号
    int t_sum=0; //当前线程完成的临时内积

    while (global_tid < n)
    {
        t_sum = a[global_tid] * b[global_tid];
        global_tid += t_n;
    }
    atomicAdd(c, t_sum);
}



int main() {
    int a[N], b[N], c;
    int *dev_a, *dev_b, *dev_c;
    for(int i=0; i<N; ++i) // 为数组a、b赋值
    {
        a[i] = i;
        b[i] = i;
    }
    cudaMalloc(&dev_a, sizeof(int) * N);
    cudaMemcpy(dev_a, a, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_b, sizeof(int) * N);
    cudaMemcpy(dev_b, b, sizeof(int) * N, cudaMemcpyHostToDevice);

    cudaMalloc(&dev_c, sizeof(int));
    cudaMemset(&dev_c, 0, sizeof(int));

    vector_dot_gpu<<<BLOCK_NUM, THREAD_NUM>>>(dev_a, dev_b, dev_c, N);

    cudaMemcpy(&c, dev_c, sizeof(int), cudaMemcpyDeviceToHost);  
    printf("c=%d\n",c);

    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);
    return 0;
}

实验：CUDA实现3*3均值滤波器

实验代码：

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

using namespace std;
using namespace std::chrono;

//相关参数设置
const int N=4096,M=4096;
const int BLOCK_NUM=32;
const int THREAD_NUM=64;
int input_img[N*M],output_img[N*M];

//初始化图像
void init_img()
{
    for(int i=0;i<N*M;i++)
        input_img[i]=i%255;
}

//cpu上实现的均值滤波器
void filter3_cpu(int w,int h,int* input_img,int *output_img)
{
    for(int i=0;i<w;i++)
        for(int j=0;j<h;j++)
        {
            int pos=i*w+j;
            if (i > 0 && (i < w - 1) && j > 0 && (j < h - 1))
            {
                int tmp=0;
                tmp+=input_img[pos];
                tmp+=input_img[pos+1];
                tmp+=input_img[pos-1];

                tmp+=input_img[pos-w];
                tmp+=input_img[pos-w+1];
                tmp+=input_img[pos-w-1];

                tmp+=input_img[pos+w];
                tmp+=input_img[pos+w+1];
                tmp+=input_img[pos+w-1];
                tmp/=9;
                output_img[pos]=tmp;
            }else{
                output_img[pos]=input_img[pos];
            }
        }
}

//gpu上实现的均值滤波器
__global__ void filter3_Gpu(int w,int h,int* input_img,int *output_img)
{
    const int tid = threadIdx.x; //当前block内的线程索引号
    const int t_n = gridDim.x*blockDim.x; // 一个grid内的线程总数
    int global_tid=blockIdx.x*blockDim.x+tid; //当前线程在整个grid中的索引号

    while(global_tid<w*h)
    {
        //根据线程id反推坐标
        int i=global_tid/h;
        int j=global_tid%h;
        int pos=i*h+j;
        //不处理边界
        if (i > 0 && (i < w - 1) && j > 0 && (j < h - 1))
        {
            int tmp=0;
            tmp+=input_img[pos];
            tmp+=input_img[pos+1];
            tmp+=input_img[pos-1];

            tmp+=input_img[pos-w];
            tmp+=input_img[pos-w+1];
            tmp+=input_img[pos-w-1];

            tmp+=input_img[pos+w];
            tmp+=input_img[pos+w+1];
            tmp+=input_img[pos+w-1];
            tmp/=9;
            output_img[pos]=tmp;
        }else{
            //边界不处理
            output_img[pos]=input_img[pos];
        }
        //线程处理下一个任务
        global_tid+=t_n;
    }
}

int main()
{
    init_img();
    auto start = system_clock::now();

    filter3_cpu(N,M,input_img,output_img);

    auto end = system_clock::now();
    auto duration = duration_cast<microseconds>(end - start);
    auto total_time=double(duration.count()) * microseconds::period::num / microseconds::period::den;
    printf("cpu耗时：%lf s\n",total_time);

    int *dev_input_img, *dev_output_img;
    start = system_clock::now();
    cudaMalloc(&dev_input_img, sizeof(int) * N*M);
    cudaMemcpy(dev_input_img, input_img, sizeof(int) * N*M, cudaMemcpyHostToDevice);
    cudaMalloc(&dev_output_img, sizeof(int) * N*M);

    filter3_Gpu<<<BLOCK_NUM, THREAD_NUM>>>(N,M,dev_input_img,dev_output_img);

    cudaMemcpy(output_img, dev_output_img, sizeof(int)*N*M, cudaMemcpyDeviceToHost);
    end = system_clock::now();
    duration = duration_cast<microseconds>(end - start);
    total_time=double(duration.count()) * microseconds::period::num / microseconds::period::den;
    printf("Gpu耗时：%lf s",total_time);
    return 0;
}