Cuda编程小作业：cudalcs

简单看完了Cuda by Example一书，其中很多内容已经落后(听别人说现在推荐还在更新的CUDA C Programming Guide)，但是也是学到了一些东西。为了巩固基础，我需要进行一些编程练习。看网上似乎没有什么简单的cuda入门编程练习，而且好多都是用来进行图形学相关的计算，我这个不会图形学也只是想体验HPC的人夹在中间略有尴尬。

联想到半年前上的NJU的OS课程Lab中包含一个将串行算法并行化的实验，虽然我的算法能力目前依旧处于一穷二白的水平，但是已经写过一次作业了，我认为还是可以比较简单地完成这项任务。在练习cuda的同时还能训练c++。

预估代码行数只有一两百行，使用cuda的情况也和使用多线程进行并行化算法有所类似，这个算法或许很适合用cuda重写。
如果不懂将LCS算法并行化的理论实现，可以去看M2的实验指南，~~我就是这么过来的，虽然我还是不太懂~~。

我们可以把 LCS 的动态规划算法 (无论是递归求解，还是迭代求解) 都看成是一个有向无环图上的计算：为了得到f(i,j)的值，我们就需要f(i - 1,j)，f(i, j -1)和f(i - 1,j - 1)的值。

我们假设主动把进行LCS的字符串长度限制为1000以内。为了简化实现。

从这里开始我们要实现的东西就和多线程编程差的远了。得益于Cuda的特性，你可以简单理解为我们创建多于str1*str2个thread来针对DP表中每个格子（暂且叫它unit）直接进行并行计算，由于几乎不用复用thread(对于需要复用thread的情况，以我的机器(3070laptop)，对本文讲述的最简单的每个block长宽固定的实现，需要两个字符串长度分别大于32*2147483647与32*65535。而在本文限制1000长度的字符串下根本不可能。而且可以进行优化，优化后不需要使用二维block，从而直接达到x维度限制)，导致最终代码实际上很简单。

准备工作

但是在开始之前，我们需要一点准备工作，由于我在VS上编程，所以会遇到一些问题，详情可以看这个VS2022进行Cuda编程中的一点坑

将代码粘贴到文件开头

#ifdef __INTELLISENSE__
#define CUDA_KERNEL_ARGS(...)
#define __CUDACC__
#else
#define CUDA_KERNEL_ARGS(...) <<< __VA_ARGS__ >>>
#endif

前置

1
2
3

constexpr auto required_block_dim = 20;//手动设定每block中每个维度有多少个thread

constexpr auto blockDIM = std::min<int>(required_block_dim, 16);

main函数行为

获取输入

unsigned int str1_len,str2_len;
std::cin >> str1_len >> str2_len;

char* str1_arr_0copy;
cudaHostAlloc( (void**)&str1_arr_0copy,sizeof(char)*str1_len, cudaHostAllocMapped);

char* str2_arr_0copy;
cudaHostAlloc((void**)&str2_arr_0copy, sizeof(char) * str2_len, cudaHostAllocMapped);

其中，我们使用cudaHostAlloc，通过cudaHostAllocMapped的flag来取得一个zero copy memory的pinned memory page。这是一种优化的手段，大概可以理解为不直接将host中的内存复制到device设备而是在device访问特定内存时，产生一个类似page fault的，再从host memory中传输数据到device。（可惜具体情况我并不知道）理论上，对于只需要单次访问的设备内存，比起在device中分配内存并从host中复制数据到device内存中，再进行访问，这种方式明显是更快速的。虽然在本例中不一定能产生性能优化：对于str1中的N个字符，每个都需要访问M次，对于str2同理。对每个str都需要访问都需要访问N*M次，所以这种实现的性能或许真的不如分配设备内存再访问。我没有考虑哪个比较强，因为我想用用新东西。不过这也导致不能直接使用cin和string来直接获取输入。

取得字符数组

1
2
3

// When reading characters, istream_iterator skips whitespace by default 
std::copy_n(std::istream_iterator<char>(std::cin), str1_len, str1_arr_0copy);
std::copy_n(std::istream_iterator<char>(std::cin), str2_len, str2_arr_0copy);

使用c++的copy_n方法从std::cin转化成的iterator读取特定个数字符存入数组。主要是简单用用c++的iterator\

在device内存中创建DP表

int* dev_DP_table;
cudaMalloc((void**)&dev_DP_table, str1_len * str2_len * sizeof(int));
//fill the table with -1 to mark the unit is not result
cudaMemset((void*)dev_DP_table, -1, str1_len * str2_len * sizeof(int));

此表不需要对cpu可见，因为当gpu计算完成后唯一有用的数值是最后一个unit的值，只需要将其复制回host内存即可。

生成设备需要使用的DP表并调用kernel函数

1 2	`dp_table dp{dev_DP_table, str1_len,str2_len , dev_str1_arr ,dev_str2_arr}; lcs_kernel CUDA_KERNEL_ARGS(blocks, threads)(dp);`

此处的CUDA_KERNEL_ARGS是前文所言文件开头define的宏，由于vs使用的intellisense无法识别核函数的尖括号调用参数会报错，我不是很喜欢，所以定义一个对intellisense来说展开后是空的宏让它不报错。

获取结果，回收内存

cudaMemcpy((void*)&output, (void*)(dp.dev_DP_table_ptr + str1_len * str2_len - 1), 1 * sizeof(int), cudaMemcpyDeviceToHost);

std::cout << output<<std::endl;
cudaFree(dev_DP_table);
cudaFreeHost(str1_arr_0copy);
cudaFreeHost(str2_arr_0copy);

像普通的memcpy那样通过移动指针来复制数据，唯一的不同是添加一个flag——cudaMemcpyDeviceToHost来从device复制数据到host。

辅助结构体实现

struct dp_table {
	volatile int* dev_DP_table_ptr;
	const int weight, height;
	const char* dev_str1;
	const char* dev_str2;

	__device__ bool update_when_available(int i,int j) {
		int skip_a = get_value(i - 1, j);
		int skip_b = get_value(i, j - 1);
		int skip_middle = get_value(i - 1, j - 1);
		//printf("[%d,%d]: %d %d %d\n", i, j, skip_a,skip_a,skip_middle);

		//no result
		if (skip_a == -1 || skip_b == -1 || skip_middle == -1) {
			return false;
		}

		int take_both = skip_middle + (dev_str1[i] == dev_str2[j]);
		//update
		set_value(i,j,max(take_both, max(skip_a, skip_b)));
		return true;
	}
	__device__ int get_value(int i, int j) {
		if(i < 0 || j < 0) return 0;
		assert(j < weight);
		assert(i < height);
		return dev_DP_table_ptr[i * weight + j];
	}
	__device__ void set_value(int i, int j, int v) {
		dev_DP_table_ptr[i * weight + j] = v;

	}
};

这是前文所述的dp_table结构体实现。让kernel函数对dp表的具体到数组的操作解耦。

这其中volatile关键字非常重要，就跟OS中多线程编程一样，在没有修饰这个指针时，指针指向的数据会被缓存，导致数据不能及时更新。在我这里，如果不加volatile，会导致单个block内threads数量的两字符串长度可正常解析，需要启动超过一个grid大小的thread就会一直卡住无法获取被修改过的数据。

kernel函数实现

__global__ void lcs_kernel(dp_table dp) {
	int j = threadIdx.x + blockDim.x * blockIdx.x;
	int i = threadIdx.y + blockDim.y * blockIdx.y;
	if (i < dp.height && j < dp.weight) {

		//spinning before updated
		while (true) {
			if (dp.update_when_available(i, j)) break;
		}

	}

	
}

通过封装辅助结构体，我尽量简化了kernel函数的实现。其中的if用于判断是否有多余的thread启动，因为我调用kernel的参数中为了确保在字符串长度小于threadsPerBlock的情况下也能启用，会启用略多于所需线程的线程数。

while语句类似于自旋锁，由于dp表中每个unit的状态（值）只会修改一次，显卡保证了基本类型的原子性读写(我暂时没有在文档中找到，只在StackOverflow中提到)，而kernel函数的线程执行过程大体是

取得某几个非自己的unit数据
数据是否已经计算出来，已经计算则向下运行，否则重新执行
将数据写入另外的属于自己的unit。

此过程中不会对一个unit产生写的race condition(因为只写一次)且不满足条件会不断自旋直到取得正确数据。因此在这里while充当了自旋锁。我们甚至不需要进行线程同步或者使用互斥锁。

通常的race condition主要由于修改数据时的环境与预期数据不同导致的，这是我的理解。而发生这种情况的原因主要是读取夹在两次写中间，产生W1-R1-W2-WR1的情况，读写了脏数据。但是由于本文中的情况是每个unit只会被修改一次，则只会有W(1)-R(1)-R(1)或者R(1)-W(2)-R(2)的情况出现，从而避免了race condition。

源码已上传至github，地址：https://github.com/dangjinghao/cuda_lcs