CUDA学习笔记（二）——核函数简单例程

先看代码，再说细节。

// kernel
__global__ void plus_array(int *x, int* y, int* z, int len){
	int index = threadIdx.x + blockIdx.x * blockDIm.x;
	if (index < len)
		z[index] = x[index] * y[index];
}

int main(){
	int num = 1 << 20;
    int buff_size = num * sizeof(int);
    // 申请host端内存
    int *x, *y, *z;
    x = (int*)malloc(buff_size);
    y = (int*)malloc(buff_size);
    z = (int*)malloc(buff_size);

    // 初始化数据
    for (int i = 0; i < num; ++i)
    {
        x[i] = 10;
        y[i] = 90;
    }

    // 申请device端内存
    int *d_x, *d_y, *d_z;
    cudaMalloc((void**)&d_x, buff_size);
    cudaMalloc((void**)&d_y, buff_size);
    cudaMalloc((void**)&d_z, buff_size);

    // 将host数据拷贝到device   cudaMemcpy(dst, src, ...)
    cudaMemcpy((void*)d_x, (void*)x, buff_size, cudaMemcpyHostToDevice);
    cudaMemcpy((void*)d_y, (void*)y, buff_size, cudaMemcpyHostToDevice);
    // 配置block,thread数量
    dim3 block_size(256);			//每个block中thread个数
    dim3 grid_size((num + block_size.x - 1) / blockSize.x);	//每个grid中，block个数
    // 执行kernel
    plus_array << < grid_size, block_size >> >(d_x, d_y, d_z, num);

    // 将device得到的结果拷贝到host 
    cudaMemcpy((void*)z, (void*)d_z, buff_size, cudaMemcpyHostToDevice);

    for (int i = 0; i < num; i += 1<<5)
        std::cout << z[i] << std::endl;

    // 释放device内存
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);
    // 释放host内存
    free(x);
    free(y);
    free(z);
    return 0;
}

cuda代码流程如下：

1. 分配host内存，并进行数据初始化；
2. 分配device内存，并从host将数据拷贝到device上；
3. 调用CUDA的核函数（kernel）在device上完成指定的运算；
4. 将device上的运算结果拷贝到host上；
5. 释放device和host上分配的内存。

核函数

上述代码中是一维数组的计算，按位求乘积，也就是两个一维矩阵叉乘。主要思想就是一个线程计算一个位置的对应元素乘积。
执行核函数之前需要先指定需要多少block, thread。

dim3 block_size(256);
dim3 grid_size((num + block_size.x - 1) / blockSize.x);

block_size是一个block中分配多少个线程，也就是一个block中，多少个线程并行计算。一般显卡最大支持1024个线程，但并不是线程越大越好。
grid_size是一个grid中申请多少个block，整除向下取证可能会存在一些数据没办法得到计算，所以会多分配一些。
这里都是按照１Ｄ的方式申请的。
执行核函数，需要在尖括号中指定好我们设置的block数，和thread数。

plus_array << < grid_size, block_size >> >(d_x, d_y, d_z);

接下来就是最重要的核函数了。

__global__ void plus_array(int *x, int* y, int* z，int len){
	int index = threadIdx.x + blockIdx.x * blockDIm.x;
	if (index < len)		// 防止内存越界
		z[index] = x[index] * y[index];
}

首先核函数在声明的时候需要加上__global__关键字。表示是在device端执行的代码。
首先需要找到当前执行这个核函数的线程，所应该操作的数据是哪一块的（会在下面内存使用的部分详细介绍），我们这里相当于一个thread只执行数组中第index个数据乘积运算。
这样我们就需要先计算出，当前的index是什么。
threadIdx, blockIdx, blockDIm是CUDA中内置的对象，由于记录当前所在位置，因为我们申请的grid, block都是1D的，所以只需要关注x的值，y,z默认都是1。
代码中grid_size是4096。计算方式如下图（偷懒引用自https://blog.youkuaiyun.com/xiaohu2022/article/details/79599947，不自己画了）

再实现一个二维矩阵的求和运算，index的计算就更为复杂了。主要想法依然是一个线程负责运算一个位置的求和。

__global__ void matrix_plus(int *src1, int *src2, int *dst, int width, int height){
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    if (x + y * width < width * height)
        dst[x + y * width] = src1[x + y * width] + src2[x + y * width];
}

void matrix_plus_2d(int *a, int *b, int width, int height) {

    cudaSetDeviceFlags(cudaDeviceMapHost);

    int *dev_src1, *dev_src2, *dev_dst, *dst, *src1, *src2;
    int byte_size = width * height * sizeof(int);
	
	// cudaHostAlloc 申请的内存是锁页内存，可以在cpu和设备端同时访问，避免相互拷贝。
    cudaHostAlloc((void**)&src1, byte_size, cudaHostAllocWriteCombined | cudaHostAllocMapped);
    cudaHostAlloc((void**)&src2, byte_size, cudaHostAllocWriteCombined | cudaHostAllocMapped);
    cudaHostAlloc((void**)&dst, byte_size, cudaHostAllocWriteCombined | cudaHostAllocMapped);
    memcpy((void*)src1, (void*)a, byte_size);
    memcpy((void*)src2, (void*)b, byte_size);
	
	// 通过引用的方式，设备端指针指向所申请的锁页内存
    cudaHostGetDevicePointer(&dev_src1, src1, 0);
    cudaHostGetDevicePointer(&dev_src2, src2, 0);
    cudaHostGetDevicePointer(&dev_dst, dst, 0);


    dim3 block_size(16, 16);
    dim3 grid_size((width / block_size.x) + 1, height / block_size.y + 1);

    matrix_plus<<<grid_size, block_size>>>(dev_src1, dev_src2, dev_dst, width, height);

	//　需要在执行完成后，同步等待所有的线程执行结束
    cudaDeviceSynchronize();

    for (int i = 0; i < width; i += 1 << 2) {
        for (int j = 0; j < height; j += 1 << 2) {

            std::cout << dst[i*width + j] << "  ";
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;

    cudaFreeHost(src1);
    cudaFreeHost(src2);
    cudaFreeHost(dst);
}

代码中使用了锁页内存的方式申请内存，具体原理等搞清楚后，在后续博客更新。基本执行顺序还是一样的。
这里需要介绍的是二维block和thread怎么访问指定位置数据，也就是计算索引的问题。

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

可见是一维基础上分别计算出x和y方向上所在的位置，x表示在x方向上左边有几个线程，y表示在y方向上方有几个线程。

再看这张结构图，block(2, 0)中thread(1, 0)，可以计算得出x=７，y=0。实际上，下标是从0开始，也就是线程在整个线程阵列中的索引坐标。
而传进去的数据是一维的，所以还需要计算，二维坐标对应一维数据的下标x + y * width，这个很简单，就不多做介绍，如果申请二维数据，这一步操作就省去了。

CMakeLists.txt

我实在linux环境下，写的代码，使用cmake进行编译，简单介绍下CMakeLists的编写。

#　主工程目录下CMakeLists
# 搜索本机CUDA环境配置
FIND_PACKAGE(CUDA REQUIRED)
if (NOT CUDA_FOUND)
    message(STATUS "CUDA not found. Project will not be built.")
endif(NOT CUDA_FOUND)
include_directories(${CUDA_INCLUDE_DIRS})

set(CMAKE_CXX_STANDARD 11)
#　cuda代码放在了cuda_code文件夹下
ADD_SUBDIRECTORY(cuda_code)

add_executable(ttest main.cpp )

target_link_libraries(ttest cuda_code_lib)


#　cuda_code的CMakeLists
message(${CMAKE_CURRENT_SOURCE_DIR})

file(GLOB cuda_code_src ${CMAKE_CURRENT_SOURCE_DIR}/*.cc )

include_directories(${CUDA_INCLUDE_DIRS})
# complie *.cu file
file(GLOB cuda_code_cu ${CMAKE_CURRENT_SOURCE_DIR}/*.cu)

# 注意：　use cuda compiler
cuda_compile(cuda_objs ${cuda_code_cu})
list(APPEND cuda_code_src ${cuda_objs} ${cuda_code_cu})

add_library(cuda_code_lib ${cuda_code_src})

target_link_libraries(cuda_code_lib ${CUDA_LIBRARIES})