《CUDA编程》3.简单CUDA程序的基本框架

本章将学习CUDA程序的基本框架，编写更加有用的CUDA程序

0 C++例子：数组相加

C++代码如下：

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h> // 包含 time.h 头文件以使用 clock()

const double EPS = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;

void add(const double* x, const double* y, double* z, const int N) {
    for (int n = 0; n < N; ++n) {
        z[n] = x[n] + y[n];
    }
}

void check(const double* z, const int N) {
    bool has_error = false;
    for (int n = 0; n < N; ++n) {
        if (fabs(z[n] - c) > EPS) {//判断浮点数是否相等时，不能用==，而应该做差，判断差值是否小于某个值
            has_error = true;
        }
    }
    printf("Has error: %d\n", has_error);
}

int main(void) {
    const int N = 100000000;
    const int M = sizeof(double) *N;

    double* x = (double*)malloc(M);
    double* y = (double*)malloc(M);
    double* z = (double*)malloc(M);

    // 记录程序开始时间
    clock_t start = clock();

    for (int n = 0; n < N; ++n) {
        x[n] = a;
        y[n] = b;
    }

    add(x, y, z, N);
    check(z, N);

    // 记录程序结束时间
    clock_t end = clock();

    // 计算并打印程序运行时间
    double cpu_time_used = ((double)(end - start)) / CLOCKS_PER_SEC;
    printf("Time used: %f seconds\n", cpu_time_used);

    free(z);
    free(x);
    free(y);

    return 0;
}

输出结果为
在这里插入图片描述
其中

    const int N = 100000000; //定义数组的长度为10的8次方
    const int M = sizeof(double) *N; // 每个数组所需的字节数

    double* x = (double*)malloc(M); // 分配内存
    double* y = (double*)malloc(M); // 分配内存
    double* z = (double*)malloc(M); // 分配内存

会创建3个长度为10的8次方的一维数组，每个数组大约占用 800 MB 的内存，总共需要约 2.4 GB 的主机内存来存储这三个数组，同样的后面的CUDA程序也需要2.4GB的GPU内存，如果不足，请自行调整

1 CUDA的基本框架

一个典型的基本框架如下图所示：
请添加图片描述
现在我们根据上述框架，将刚刚的C++程序编写为CUDA程序：

#include <cuda.h>
#include <cuda_runtime.h>
#include <math.h>
#include <stdio.h>
#include <time.h> // 包含 time.h 头文件以使用 clock()

const double EPS = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;

// 希望 add 函数在 GPU 上执行
__global__ void add(const double *x, const double *y, double *z);
void check(const double* z, const int N);

int main(void) {
    const int N = 100000000; // 定义数组的长度为 10 的 8 次方
    const int M = sizeof(double) * N; // 每个数组所需的字节数

    // 分配host内存
    double* h_x = (double*)malloc(M); 
    double* h_y = (double*)malloc(M); 
    double* h_z = (double*)malloc(M); 

    // 记录程序开始时间
    clock_t start = clock();

    for (int n = 0; n < N; ++n) {
        h_x[n] = a;
        h_y[n] = b;
    }
	
	//分配device内存
    double* d_x, * d_y, * d_z;
    cudaMalloc((void**)&d_x, M);
    cudaMalloc((void**)&d_y, M);
    cudaMalloc((void**)&d_z, M);

    // 将数据从主机复制到设备上
    cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice);
    cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice);

    const int block_size = 128;
    // 计算网格尺寸，确保所有元素都能被处理
    const int grid_size = (N + block_size - 1) / block_size; 

    // 调用内核函数在设备中进行计算
    add <<<grid_size, block_size>> > (d_x, d_y, d_z);

    // 将计算结果从设备复制回主机
    cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost);
    check(h_z, N);

    // 记录程序结束时间
    clock_t end = clock();

    // 计算并打印程序运行时间
    double cpu_time_used = ((double)(end - start)) / CLOCKS_PER_SEC;
    printf("Time used: %f seconds\n", cpu_time_used);
	
	// 释放内存
    free(h_x);
    free(h_y);
    free(h_z);
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);
    return 0;
}

__global__ void add(const double* x, const double* y, double* z) {
    const int n = blockIdx.x * blockDim.x + threadIdx.x;
    z[n] = x[n] + y[n];
}

void check(const double* z, const int N) {
    bool has_error = false;
    for (int n = 0; n < N; ++n) {
        if (fabs(z[n] - c) > EPS) {
            has_error = true;
        }
    }
    printf("Has error: %d\n", has_error);
}

输出结果如下：
请添加图片描述

1.1 解释grid_size的计算

const int block_size = 128;
// 计算网格尺寸，确保所有元素都能被处理
const int grid_size = (N + block_size - 1) / block_size;

其中

N：要处理的元素总数
block_size：每个block中的线程数，这里是128，可以根据任务进行调整，最大是1024

公式(N + block_size - 1) / block_size是确保N不能被block_size 整除，即确保有足够的线程块来处理元素，下面举例说明：

整除情况：例如N=1024，block_size=128，则1024/128=8，所以需要8个线程块，则grid_size设置为8
不能整除：例如N=1025，block_size=128，则1025/128=8余1，所以还需要多一个线程块来处理，但是grid_size是int类型，所以此时也是取8。
为了避免这种情况，通过加 block_size - 1 来确保 N 在除以 block_size 时总是向上取整。

1.2 隐形的device初始化

在CUAD的runtime的API中，没有明显的初始化设备的函数，因为在第一次调用一个设和设备的管理及版本查询功能无关的runtime的API时，device将自动初始化

1.3 device内存的分配与释放

在代码中，我们在hsot和device中分别定义了3个数组，也分别分配了内存和显存
在这里插入图片描述
大家发现第一个指针是一个双重指针，是因为cudaMalloc()的函数原型是

第一个参数address是待分配设备内存的指针，但内存本身就是一个指针，所以该参数其实是指针的指针，所以使用双重指针
size是待分配内存的字节数
返回值是一个错误代号。如果调用成功，返回cudaSuccess，否则返回一个错误代码。

看下面这段代码
在这里插入图片描述
虽然double类型占用是8字节，但是为了程序的健壮性和可移植性（有的设备可能不是），所以这里还是使用sizeof(double)

1.4 host和device之间数据的传递

这里使用的是cudaMemcpy()函数，原型是：
在这里插入图片描述

dst：目标地址
src：源地址
count：复制数据的字节数
kind：标志，取以下几个值
–cudaMemcpyHostToHost，表示从主机复制到主机
–cudaMemcpyHostToDevice，表示从主机复制到设备
–cudaMemcpyDeviceToHost，表示从设备复制到主机
–cudaMemcpyDeviceToDevice，表示从设备复制到设备
–cudaMemcpyDefault，表示根据dst和src所指的地址自动判断，要求系统有统一寻址功能（unified virtual addressing）
返回值是一个错误代号，如果调用成功，返回cudaSuccess
该函数将一定字节数的数据从源地址所指的缓冲区复制到目标地址所指的缓冲区

所以下列代码的意义是：

cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice);
//将h_x指向的主机内存中的数据，复制到d_x指向的设备内存中去

1.5 核函数中数据与线程的对应

观察下列代码
在这里插入图片描述
使用了一个一维线程块，一共有【10的8次方/128】个线程块，每个线程块有128个线程

再观察下列函数：

这个是C++中写的函数
在这里插入图片描述
这个是CUDA程序中写在device中的函数

可以发现，将主机函数修改为设备中的核函数，只需要去掉循环即可。

因为在主机函数中，需要依次对数组的元素进行操作，所以会写一个循环。
但是在核函数中，我们使用“单指令— —多线程”的方式进行编写，即每个元素的操作都由单独的一个线程来计算，所以不需要循环，会进行并行执行

注意： 也可以使用循环，只需要设置为<<<1, 1>>>，就可以用一个线程去调用核函数计算

2 核函数的要求

2.1 编写核函数时要注意：

核函数的返回类型必须是void
必须使用__global__，也可以加上其他C++的限定符，如static，次序任意
函数名无特殊要求，支持重载
不支持可变数量的参数列表，参数的个数必须确定
核函数不能成为一个类的成员，通常是用一个包装函数调用核函数，然后将包装函数定义为类的成员
在计算能力3.5之前的，核函数之间不能相互调用
可以向核函数传递非指针变量（如int N），其内容每个线程可见（可使用可访问）
除非使用统一内存编程机制（将在第12章介绍），否则传给核函数的数组（指针）必须指向设备内存

2.2 核函数中if语句的必要性

在这里插入图片描述
该代码中没有使用if语句，但由于线程数是比元素总数多的，所以应该添加一个if语句，来防止数组越界，并且节约不必要的线程计算，修改代码如下：

运行结果如下，比不加if语句之前快，因为多余的线程没有执行计算操作
在这里插入图片描述

3 自定义设备函数

核函数可以调用不带执行配置的自定义函数，这样的自定义函数称为设备函数（device function），它是在设备中执行，并在设备中被调用的。

与之相比，核函数是在设备中执行，但在主机端被调用的。

3.1 标识符

用__global__修饰的函数称为核函数，一般由主机调用，在设备中执行。如果使用动态并行，则也可以在核函数中调用自己或其他核函数。
用__device__修饰的函数叫称为设备函数，只能被核函数或其他设备函数调用，在设备中执行。
用__host__修饰的函数就是主机端的普通C++函数，在主机中被调用，在主机中执行。对于主机端的函数，该修饰符可省略。因为有时可以用__host__和__device__同时修饰一个函数，使得该函数既是一个C++中的普通函数，又是一个设备函数。这样做可以减少冗余代码。编译器将针对主机和设备分别编译该函数
不能同时用__device__和__global__修饰一个函数，即不能将一个函数同时定义为设备函数和核函数。
不能同时用__host__和__global__修饰一个函数，即不能将一个函数同时定义为主机函数和核函数。
编译器决定把设备函数当作内联函数（inline function）或非内联函数，但可以用修饰符__noinline__建议一个设备函数为非内联函数（编译器不一定接受），也可以用修饰符__forceinline__建议一个设备函数为内联函数

3.2 例子：为数组相加的核函数定义一个设备函数

（1）返回值

#include <cuda.h>
#include <cuda_runtime.h>
#include <math.h>
#include <stdio.h>
#include <time.h> // 包含 time.h 头文件以使用 clock()

const double EPS = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;

// 希望 add 函数在 GPU 上执行
__global__ void add(const double *x, const double *y, double *z,int N);
__device__ double add1_device(const double x, const double y);
void check(const double* z, const int N);

int main(void) {
    const int N = 100000000; // 定义数组的长度为 10 的 8 次方
    const int M = sizeof(double) * N; // 每个数组所需的字节数

    // 在 host 中分配内存
    double* h_x = (double*)malloc(M); // 分配内存
    double* h_y = (double*)malloc(M); // 分配内存
    double* h_z = (double*)malloc(M); // 分配内存

    // 记录程序开始时间
    clock_t start = clock();

    for (int n = 0; n < N; ++n) {
        h_x[n] = a;
        h_y[n] = b;
    }

    double* d_x, * d_y, * d_z;
    cudaMalloc((void**)&d_x, M);
    cudaMalloc((void**)&d_y, M);
    cudaMalloc((void**)&d_z, M);

    // 将主机上的数据复制到设备上
    cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice);
    cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice);

    const int block_size = 128;
    // 计算网格尺寸，确保所有元素都能被处理
    const int grid_size = (N + block_size - 1) / block_size; 

    // 调用内核函数来执行加法
    add <<<grid_size, block_size>> > (d_x, d_y, d_z,N);

    // 将计算结果从设备复制回主机
    cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost);
    check(h_z, N);

    // 记录程序结束时间
    clock_t end = clock();

    // 计算并打印程序运行时间
    double cpu_time_used = ((double)(end - start)) / CLOCKS_PER_SEC;
    printf("Time used: %f seconds\n", cpu_time_used);

    free(h_x);
    free(h_y);
    free(h_z);
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);

    return 0;
}

//返回值
__device__ double add1_device(const double x, const double y) {
    return x + y;
}

__global__ void add(const double* x, const double* y, double* z,int N) {
    const int n = blockIdx.x * blockDim.x + threadIdx.x;//会飘红，不影响运行
    // 添加边界检查，确保索引不越界
    if (n < N) {
        z[n] = add1_device(x[n], y[n]);
    }
}

void check(const double* z, const int N) {
    bool has_error = false;
    for (int n = 0; n < N; ++n) {
        if (fabs(z[n] - c) > EPS) {
            has_error = true;
        }
    }
    printf("Has error: %d\n", has_error);
}

（2）指针

#include <cuda.h>
#include <cuda_runtime.h>
#include <math.h>
#include <stdio.h>
#include <time.h> // 包含 time.h 头文件以使用 clock()

const double EPS = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;

// 希望 add 函数在 GPU 上执行
__global__ void add(const double *x, const double *y, double *z,int N);
__device__ void add2_device(const double x, const double y,double *z);
void check(const double* z, const int N);

int main(void) {
    const int N = 100000000; // 定义数组的长度为 10 的 8 次方
    const int M = sizeof(double) * N; // 每个数组所需的字节数

    // 在 host 中分配内存
    double* h_x = (double*)malloc(M); // 分配内存
    double* h_y = (double*)malloc(M); // 分配内存
    double* h_z = (double*)malloc(M); // 分配内存

    // 记录程序开始时间
    clock_t start = clock();

    for (int n = 0; n < N; ++n) {
        h_x[n] = a;
        h_y[n] = b;
    }

    double* d_x, * d_y, * d_z;
    cudaMalloc((void**)&d_x, M);
    cudaMalloc((void**)&d_y, M);
    cudaMalloc((void**)&d_z, M);

    // 将主机上的数据复制到设备上
    cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice);
    cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice);

    const int block_size = 128;
    // 计算网格尺寸，确保所有元素都能被处理
    const int grid_size = (N + block_size - 1) / block_size; 

    // 调用内核函数来执行加法
    add <<<grid_size, block_size>> > (d_x, d_y, d_z,N);

    // 将计算结果从设备复制回主机
    cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost);
    check(h_z, N);

    // 记录程序结束时间
    clock_t end = clock();

    // 计算并打印程序运行时间
    double cpu_time_used = ((double)(end - start)) / CLOCKS_PER_SEC;
    printf("Time used: %f seconds\n", cpu_time_used);

    free(h_x);
    free(h_y);
    free(h_z);
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);

    return 0;
}

//指针
__device__ void add2_device(const double x, const double y, double* z) {
    *z=x+y;
}

__global__ void add(const double* x, const double* y, double* z,int N) {
    const int n = blockIdx.x * blockDim.x + threadIdx.x;//会飘红，不影响运行
    // 添加边界检查，确保索引不越界
    if (n < N) {
        add2_device(x[n], y[n], &z[n]);
    }
}

void check(const double* z, const int N) {
    bool has_error = false;
    for (int n = 0; n < N; ++n) {
        if (fabs(z[n] - c) > EPS) {
            has_error = true;
        }
    }
    printf("Has error: %d\n", has_error);
}

（3）引用

#include <cuda.h>
#include <cuda_runtime.h>
#include <math.h>
#include <stdio.h>
#include <time.h> // 包含 time.h 头文件以使用 clock()

const double EPS = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;

// 希望 add 函数在 GPU 上执行
__global__ void add(const double *x, const double *y, double *z,int N);
__device__ void add3_device(const double x, const double y,double &z);
void check(const double* z, const int N);

int main(void) {
    const int N = 100000000; // 定义数组的长度为 10 的 8 次方
    const int M = sizeof(double) * N; // 每个数组所需的字节数

    // 在 host 中分配内存
    double* h_x = (double*)malloc(M); // 分配内存
    double* h_y = (double*)malloc(M); // 分配内存
    double* h_z = (double*)malloc(M); // 分配内存

    // 记录程序开始时间
    clock_t start = clock();

    for (int n = 0; n < N; ++n) {
        h_x[n] = a;
        h_y[n] = b;
    }

    double* d_x, * d_y, * d_z;
    cudaMalloc((void**)&d_x, M);
    cudaMalloc((void**)&d_y, M);
    cudaMalloc((void**)&d_z, M);

    // 将主机上的数据复制到设备上
    cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice);
    cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice);

    const int block_size = 128;
    // 计算网格尺寸，确保所有元素都能被处理
    const int grid_size = (N + block_size - 1) / block_size; 

    // 调用内核函数来执行加法
    add <<<grid_size, block_size>> > (d_x, d_y, d_z,N);

    // 将计算结果从设备复制回主机
    cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost);
    check(h_z, N);

    // 记录程序结束时间
    clock_t end = clock();

    // 计算并打印程序运行时间
    double cpu_time_used = ((double)(end - start)) / CLOCKS_PER_SEC;
    printf("Time used: %f seconds\n", cpu_time_used);

    free(h_x);
    free(h_y);
    free(h_z);
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);

    return 0;
}

//指针
__device__ void add3_device(const double x, const double y, double& z) {
    z=x+y;
}

__global__ void add(const double* x, const double* y, double* z,int N) {
    const int n = blockIdx.x * blockDim.x + threadIdx.x;//会飘红，不影响运行
    // 添加边界检查，确保索引不越界
    if (n < N) {
        add3_device(x[n], y[n], z[n]);
    }
}

void check(const double* z, const int N) {
    bool has_error = false;
    for (int n = 0; n < N; ++n) {
        if (fabs(z[n] - c) > EPS) {
            has_error = true;
        }
    }
    printf("Has error: %d\n", has_error);
}