基于CUDA的GPU并行计算技术实现网课课表编排

这篇文章是用来填这个坑的：https://blog.csdn.net/xinew4712/article/details/108276264

上篇文末设想的是用天灾和定向改造机制来提高排课运算的效率，结果并不尽如人意。虽然如此，我还是会把天灾和基因改造算法帖出来，抛砖引玉吧。

而为了提高效率，最终采用的是基于CUDA的GPU并行计算技术，这也是个不小的坑，而且不论中文环境还是github都没有类似的样本先例，github上有两个基于CUDA的时间表算法代码，使用的是假设的损失函数做核函数，没有实际应用价值。

可能时间表问题本身并不适合做GPU并行计算，该问题形成的矩阵是变长的，而并行算法擅长的图像处理，都是使用定长（指定分辨率）的图片。

而本文还是用基于CUDA的GPU并行技术把网课课表编排实现了，水平有限，其中有两个算法写好了核函数，可是长时间运行，总是报错，并行代码的bug太难查了，索性采用串并结合的方案把任务跑下来了。

这是可以解决实际问题的代码，已经在实际教学中用于课表编排，文末会提供源码和真实数据（已脱敏），也请各位大牛批评指点，隐藏其中的bug是怎么回事？

即便是串并结合，效率也大大高于纯CPU代码，还是有一定的应用价值的。 

先把上篇文章的坑填了，天灾机制：

            # 在取得精英种群之后，50代以内3次不进化，启动天灾，消灭精英，形成突变，为种群发展提供更多可能。
            # 天灾机制有待研究，30%的精英率，经过4代的繁衍，就可以让精英的后代布满整个种群，
            # 百代之后的天灾消灭不了精英，种群发展会进入死胡同。应该全部杀掉，重新初始化。

            if i < 50 and noEvolution > 2:
                noEvolution = 1
                print('进程号：{} | 迭代: {} | 发生天灾'.format(os.getpid(), i + 1))
                newPopulation = []
                psize = [j for j in range(0, self.popsize)]
                #list求差，找出所有非精英个体，然后取序号前30。
                c = list(set(psize).difference(set(eliteIndex)))
                for k in range(0, self.elite):
                    newPopulation.append(self.population[c[k]])
            if i > 300 and noEvolution > 100 and bestScore > 30:
                print('进程号：{} | 迭代: {} | 重新进化'.format(os.getpid(), i + 1))
                # 开启新的进化。
                break

在不改变损失值的情况下，生成所有可能的方案，填充种群，让这样定向改造后的种群，自由繁衍：

def test(schedules):
    #对每一个元素执行一次，生成所有可以定向改造的方向。
    index1 = schedulesId.index(input_str)
    print("原时段：" + schedules[index1].liveTime)
    for i in range(len(schedulesTime)):
        conflict = 0
        schedules[index1].liveTime = schedulesTime[i]
        for m in range(0, len(schedules) - 1):
            for n in range(m + 1, len(schedules)):                
                if schedules[m].teacherId == schedules[n].teacherId and str(schedules[m].liveTime) == str(schedules[n].liveTime):
                    conflict += 1
                if str(schedules[m].liveTime) == str(schedules[n].liveTime):
                    conflict += len(schedules[m].classId) + len(schedules[n].classId) - len(list(set(schedules[m].classId +
                                                                                                     schedules[n].classId)))
        print("时段"+str(i+1)+":"+str(schedulesTime[i])+"，冲突数:"+str(conflict))

下面开始本文的正文。

 

• 一、问题的引出。

为什么要考虑用GPU并行计算来解决排课问题，主要原因是CPU慢。

第一次用遗传算法解决排课问题时，补考报名还没有开放，所有上网课的学生都是正考的，为正考学生排课，有简便的解决思路：即用班级号代替学号检查冲突，因为同一班级的所有学生选课是一样的，只要保证班号不存在课程冲突，学号也一定不会冲突，而班号只有1万多个，课程200来门，整体计算量并不大，传统的CPU串行编程也能在比较短的时间内算出可行解来。

而当补考的学生也开始上网课，不少学生开始报怨课程冲突，后面再次排课时，有必要考虑补考学生的正常学习需要。当然，混学分的不在考虑范围内，有的学生一个学期补考10来门，为了计算简便是不需要考虑的，设计目标是要求补考科目在4门以下的学生，保证其课表无冲突。这时，一个自然班每个学生选课开始出现差异，用班号来检查冲突不可行，唯有靠学号了。

统计一下选课人数，近200门课，近50万条选课记录，运算复杂度要爆炸了。虽然理论上，使用一样的算法，CPU也一定能算出可行解，但时间的消耗十分巨大，后面会有对比图。

我们知道，一般的电脑，CPU有十几二十个核心算是比较不错的了。而一般的GPU，有一两千的核心都是很普通很平常的，虽然GPU核心没有CPU核心那么强大，但它胜在数量多，这里举个形象的例子：20个数学老师和2000个小学生PK，同时运算2万道100以内加减法，请问，谁先做完？答案不言而喻吧。这就是考虑使用GPU进行科学计算的理论基础了。而程序猿的任务就是把复杂问题简单化，把具体任务切割成并行计算的一个一个小任务，同时丢给几千个GPU核心去计算，再收集处理返回结果。大多数情况下，传统编程中的for循环都适合改造成并行代码，比如第3次循环的结果并不影响第5次的运算，这样的场景都适合切割任务后交给GPU并行运算。

• 二、初识并行计算。

对于编写串行代码20年的老程序猿，串行思想根深蒂固，转变思想是最难的，要跳出习以为常的舒适思维框架，重新以并行的角度思考问题。想简单加个装饰器就并行化，那样的想法太幼稚了，代码基本是要重写的，由于经验不足，重写了还不止一次。

其实，在去年研究计算机深度学习时，就接触过GPU并行计算，但那基本上是使用封装好的框架，是基于Tensorflow的并行计算，自己并不清楚底层如何运作，数据是怎么在CPU与GPU之间流动，全局显存、共享显存有什么区别，如何基于线程、块、网格分割数据使代码更高效运行，自己的显卡有多少流处理器，有多少cuda核，每个块有多少缓存，一概不知。这次使用numba库基于cuda开发小程序，算是把这些弄清楚了，甚至还写了一些C代码进行性能对比。

• 三、代码重温。

GPU端运行的代码，称为核函数。核函数是比较娇气的，有两样东西她是不支持的：这也不支持，那也不支持。

1、自定义对象是用不了的，甚至连list都不行，所有运算都基于numpy矩阵，那就直接从csv中读取数据：

# 读取课程教师关系数据,有表头，去掉。
courseTeacher = np.loadtxt("222.csv", delimiter=",", skiprows=1, dtype=np.int32)
# 读取课程与班级（学号）对应关系数据。注意生成数据时取classId的“id”。对应关系为班级时，只考虑正考。对应关系为学号时，考虑了补考。
courseStudent = np.loadtxt("555.csv", delimiter=",", dtype=np.int32)

2、由于最简单的排序都不支持，在生成数据源时，在excel中顺手给它排个序，程序处理时要高效方便的多。

3、随机函数不支持。我也不知道怎么在核函数中写生成随机数算法（据说可以），那就在CPU中生成好，传给她。这时要培养一个新习惯，CPU端变量加上“_host”，GPU端变量加上”_device”，不然容易搞混了。

randTwinsArray_host = np.random.randint(0, Maxrow, size=[n, 2])
randTwinsArray_device = cuda.to_device(randTwinsArray_host)
randIntArray_host = np.random.randint(0, Esize, size=[n, ])
randIntArray_device = cuda.to_device(randIntArray_host)

4、变量定义支持，但要非常小心。

 

Local Memory、Shared Memory分别是多少，要根据自己的显卡跑代码查询一下，矩阵大了，不小心就溢出了。

def query_device():
    drv.init()
    print('CUDA device query (PyCUDA version) \n')
    print(f'Detected {drv.Device.count()} CUDA Capable device(s) \n')
    for i in range(drv.Device.count()):

        gpu_device = drv.Device(i)
        print(f'Device {i}: {gpu_device.name()}')
        compute_capability = float( '%d.%d' % gpu_device.compute_capability() )
        print(f'\t Compute Capability: {compute_capability}')
        print(f'\t Total Memory: {gpu_device.total_memory()//(1024**2)} megabytes')

        # The following will give us all remaining device attributes as seen
        # in the original deviceQuery.
        # We set up a dictionary as such so that we can easily index
        # the values using a string descriptor.

        device_attributes_tuples = gpu_device.get_attributes().items()
        device_attributes = {}

        for k, v in device_attributes_tuples:
            device_attributes[str(k)] = v

        num_mp = device_attributes['MULTIPROCESSOR_COUNT']

        # Cores per multiprocessor is not reported by the GPU!
        # We must use a lookup table based on compute capability.
        # See the following:
        # http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#compute-capabilities

        cuda_cores_per_mp = { 7.5 : 64, 5.1 : 128, 5.2 : 128, 6.0 : 64, 6.1 : 128, 6.2 : 128}[compute_capability]

        print(f'\t ({num_mp}) Multiprocessors, ({cuda_cores_per_mp}) CUDA Cores / Multiprocessor: {num_mp*cuda_cores_per_mp} CUDA Cores')

        device_attributes.pop('MULTIPROCESSOR_COUNT')

        for k in device_attributes.keys():
            print(f'\t {k}: {device_attributes[k]}')

5、CPU与GPU是异步计算的，CPU把任务交待下去了，GPU就老老实实算，算没算完，CPU不管的，要想等待GPU的结果再做下一步打算，要同步一下。

# 等待所有并发线程计算结束
cuda.synchronize()

6、越是高频、简单的运算，在GPU端并发执行，效果越好，在核函数中跑一些小小的循环，可以高效复用线程，毕竟线程的创建和销毁也要消耗资源的，下面是在核函数中跑损失函数：

@cuda.jit
def gpu_cost(cs, ct, conflicts, n):
    pos = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x
    # if pos == 0 :
    #     sA = cuda.shared.array(shape=(175, 3), dtype=np.int32)
    #     sA = ct
    # cuda.syncthreads()
    if pos < n:
        conflicts[pos] = 0
        # 并行计算，用该线程索引所在元素，与该元素后面的每一个元素做对比
        for j in range(pos + 1, n):
            # 学号相同则判断是否冲突，不同不冲突
            if cs[pos][1] == cs[j][1]:
                # 历遍ct表，找courseId对应的时间
                for i in ct:
                    # 该线程索引所在元素的courseId与ct相同，则改courseId为LiveTimeId
                    if cs[pos][0] == i[0]:
                        cs[pos][0] = i[2] * 100000000 + cs[pos][1]
                    # 对后面的每一个元素转换courseId到LiveTimeId
                    if cs[j][0] == i[0]:
                        cs[j][0] = i[2] * 100000000 + cs[j][1]
                 # 如果有相等，那就是冲突了，返回该线程值为1
                if cs[pos][0] == cs[j][0]:
                    conflicts[pos] = 1
                    break
            else:
                break

7、同样的，交叉和变异函数也是高频调用的，也写好了核函数，它能运行，但经不住时间考验，跑着跑着就unknow ERROR了，还没查出原因：

@cuda.jit
def gpu_cross(newPopulation_device,
             resultCrossPopulation_device,
             Maxrow,
             randTwinsArray_device,
             randIntArray_device,
             n_cross):

    pos, row = cuda.grid(2)
    if pos < n_cross and row < Maxrow:
        new_row = randIntArray_device[pos] * Maxrow + row
        resultCrossPopulation_device[pos][row][0] = newPopulation_device[new_row][0]
        resultCrossPopulation_device[pos][row][1] = newPopulation_device[new_row][1]
        resultCrossPopulation_device[pos][row][2] = newPopulation_device[new_row][2]
        cuda.syncthreads()
        if row == 0:
            index_a = randIntArray_device[pos] * Maxrow + randTwinsArray_device[pos][0]
            index_b = randIntArray_device[pos] * Maxrow + randTwinsArray_device[pos][1]
            tmp = resultCrossPopulation_device[pos][index_a][2]
            resultCrossPopulation_device[pos][index_a][2] = resultCrossPopulation_device[pos][index_b][2]
            resultCrossPopulation_device[pos][index_b][2] = tmp
@cuda.jit
def gpu_mutate(newPopulation_device,
             resultMutprobPopulation_device,
             liveTime,
             Maxrow,
             randIntArray_device,
             randRandArray_device,
             n_mutprob):
    pos, row = cuda.grid(2)
    if pos < n_mutprob and row < Maxrow:
        new_row = randIntArray_device[pos] * Maxrow + row
        resultMutprobPopulation_device[pos][row][0] = newPopulation_device[new_row][0]
        resultMutprobPopulation_device[pos][row][1] = newPopulation_device[new_row][1]
        resultMutprobPopulation_device[pos][row][2] = newPopulation_device[new_row][2]
        cuda.syncthreads()
        mut = resultMutprobPopulation_device[pos][row][2]
        if randRandArray_device[pos] > 0.5:
            if mut < liveTime:
                mut += 1
            else:
                mut -= 1
        else:
            if mut - 1 > 0:
                mut -= 1
            else:
                mut += 1
        resultMutprobPopulation_device[pos][row][2] = mut
        cuda.syncthreads()

8、流水线控制。CPU与GPU运算是异步的，GPU与总线传输也是异步的，GPU没有必要等待所有数据传输完再开始运算，数据传输一部分就可以开始运算了，这就需要流水线控制，把一个具体任务分配到不同的流，默认流与多流对比：

代码实现：

# 创建x个cuda stream
stream_list = list()
for i in range(0, number_of_streams):
    stream = cuda.stream()
    stream_list.append(stream)
# 定义每个块的线程数
threads_per_block = 64
# 每个stream的处理的数据变为原来的1/x
blocks_per_grid = math.ceil(segment_size / threads_per_block)
streams_gpu_result = np.empty(n).astype(np.int32)
streams_out_device = cuda.device_array(segment_size)

# 启动多个stream
for i in range(0, number_of_streams):
    # 传入不同的参数，让函数在不同的流执行
    # 将数据复制到显存
    x_i_device = cuda.to_device(cs[i * segment_size: (i + 1) * segment_size], stream=stream_list[i])
    y_i_device = cuda.to_device(ct, stream=stream_list[i])
    gpu_cost[blocks_per_grid, threads_per_block, stream_list[i]](
        x_i_device,
        y_i_device,
        streams_out_device,
        segment_size)
    # 将结果从显存复制回内存
    streams_out_device.copy_to_host(stream=stream_list[i])
    streams_gpu_result[i * segment_size: (i + 1) * segment_size] = streams_out_device.copy_to_host(stream=stream_list[i])
# 等待所有并发线程计算结束
cuda.synchronize()

9、变长矩阵的处理。定义好流的数量，就要分割数据，而切割时很可能把一个学号切到2个流中，这会导致冲突漏算，要手工填充无意义数据，保证同一学号完整地在同一流中，而总体上，还要把长尾进行特殊处理，这点规模的小数据，在CPU端算一下就好了，然后把结果与GPU返回的合并计算。

# 定义list，存储每个个体的冲突值
conflicts = []
# 考一份cs出来备用，避免修饰过程对原数据造成损坏。
cs = copy.deepcopy(courseStudent)
# 取得cs的原始长度
n = len(cs)
# 准备使用多少个stream，分流不能用于直接求冲突数。数据截断后会造成漏算。
number_of_streams = 5
# 每个stream处理的数据量为原来的 1/x
# 符号//得到一个整数结果，不能整除的部分不流入GPU并行，在CPU端计算
segment_size = n // number_of_streams
# 初始segment_size定了就不能再根据cs长度变化而变化了。
for k in range(0, number_of_streams):
    # print(cs[segment_size*(k+1)])
    # print("CS:" + str(len(cs))+","+str(k*segment_size)+"-"+str((k+1)*segment_size))
    # 判断断点，某流的最后一个元素与下一流的第一个元素是否相同，不同说明运气不错，相同就需要特别处理。
    # 有调试的需要，将流设为1时，本段代码要注释，越界
    if cs[segment_size * (k+1)][1] == cs[segment_size * (k+1)+1][1]:
        i = 1
        # 只要流尾部有相同数据，就循环向前继续找，直到找到不相同的那个，算出index，插入虚拟数据，将真实学号完整挤入下一流。
        while cs[segment_size * (k+1) - i][1] == cs[segment_size * (k+1)][1]:
            i += 1
        for j in range(0, i):
            # 修饰数据，插入i个不存在的学号和课程号，确保不出现某个学号被切到2个流的情况。不存在的课程号没有匹配时段，不造成冲突。
            cs = np.insert(cs, [segment_size * (k+1) + 1 - i], [888000+j, 11111111], axis=0)

# 生成去重码，将时段码乘以100000000，加上教师学生id，去重码相同，表示某时段，有教师或学生课程冲突。
clip_ct[:, 0] = clip_ct[:, 2] * 100000000 + clip_ct[:, 1]
# cs[:, 0] = cs[:, 2] * 100000000 + cs[:, 1]  转用GPU并发处理

# 切片，变一维
clip_ct = clip_ct[:, 0]
# cs = cs[:, 0]
# cs = cs.flatten().astype(np.int32)
# 0.002  cs部分转用GPU并发处理

# 去重
# uniques = np.unique(clip_ct, axis=0)
conflict_ct = len(clip_ct) - quick_remvoal(clip_ct)

# 最后，处理尾巴数据,那些不能整除的部分,尾巴一定会有，但不会太长，适合在cpu端处理
tail = cpu_cost(cs[number_of_streams * segment_size:len(cs)-number_of_streams * segment_size], ct)
# print("GPU 结果：" + str(sum(streams_gpu_result)))
conflicts.append(sum(streams_gpu_result)+tail+conflict_ct)

• 四、运算效率对比及小结

正式运算时可以用好断点机制，用一个略小于理论值的时段数配合高变异率跑上几百轮，再从断点中选优，慢慢增加时段数，配合低变异率可以快速接近最终解。

同时，业务上也有可优化的地方，比如一些选课人数过万的大课，把这样的课剔除出来，人工排到一个时段，程序运算的复杂度将呈指数下降，甚至几个小时就能算出多套可行方案来。

其它的似乎没什么好说的了。并行代码的bug很难查，唯有尽量少写bug，统一构思，多动笔，画好流图再coding，花这么大精力把代码并行化，实际效果怎么样呢？

我以1万轮运算为标准，毕竟该算法是靠迭代次数的堆积来逼近最终解的。

同样的算法在工作站上调用GPU跑：约10.5小时跑完，CPU/GPU负荷如图：

在服务器上用CPU跑实在太慢了，最终没有等它跑完：动用6台服务器共60核CPU，数据分块跑，花了88小时跑完才跑完2500轮，折合1万轮要350小时。CPU负荷如图：

 

服务器的CPU比较老了，跟工作站的GPU对比似乎不太公平，那么在工作站上用CPU跑再试一下：约76小时跑完。确实比服务器要快得多了。

当然，代码还有很多有待优化改进的地方，比如Shared Memory没有好好利用，部分算法没有完全并行化，有BUG待排查，突击实现功能导致代码写得比较乱。即便这样，整体上有8、9倍的速度提升，还是比较满意了。

附上源码，包含真实选课和教学安排数据，还赠送断点数据，部分运行日志。

https://download.csdn.net/download/xinew4712/13183427