如果这篇博客帮助到你，可以请我喝一杯咖啡~
CC BY 4.0 （除特别声明或转载文章外）
用 MPI 和 CUDA 实现链表求序（linked list ranking）算法，并与串行程序比较。
- 输入：一个 N 元素的单向链表，该链表放在一个长度为 N 的向量中。链表的每一个元素（除链尾的一个）都有一个后继的数组下标。链尾的元素的后继为
-1。 - 输出：定义链表元素的序号（rank）为「其后继的总数」，例如链头的序号为 N-1、链尾的序号为 0。向量中的元素不一定按照链表的顺序存放。要求算出每一个向量元素的序号。
实验环境
实验在老师提供的计算集群的一个节点上进行。单节点的显卡配置如下：
$ nvdia-smi
Mon Dec 2 08:38:49 2019
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 410.48 Driver Version: 410.48 |
|-------------------------------+----------------------+----------------------+
| GPU Name Persistence-M| Bus-Id Disp.A | Volatile Uncorr. ECC |
| Fan Temp Perf Pwr:Usage/Cap| Memory-Usage | GPU-Util Compute M. |
|===============================+======================+======================|
| 0 Tesla V100-PCIE... On | 00000000:3B:00.0 Off | 0 |
| N/A 30C P0 24W / 250W | 0MiB / 16130MiB | 0% Default |
+-------------------------------+----------------------+----------------------+
+-----------------------------------------------------------------------------+
| Processes: GPU Memory |
| GPU PID Type Process name Usage |
|=============================================================================|
| No running processes found |
+-----------------------------------------------------------------------------+
实验过程与分析
List Ranking 是一个比较经典的研究并行算法设计的问题，因为这个问题的并行算法相对于串行版本很难有十分明显的速度改进。但是并行化的好处并不仅仅是提速。
比如，单机上的内存总量很多情况下是不足够的，而多机并行的时候可以使用更大的内存。针对这一问题，并行化的优点是可以通过扩充硬件的规模从而扩充能够处理问题的规模，这是串行程序比不了的。
调度脚本list_ranking.pbs
这里由于临近期末，集群资源比较紧张，所有的调度都只占用集群上的一个节点进行。MPI 使用一台主机上的 32 个核进行并行；CUDA 使用节点上的单张 v100 显卡。为方便起见，我将所有的测试代码写进了一个测试脚本list_ranking.pbs中，这样直接对这个脚本进行调度运行即可自动运行完整实验并进行计时。
#PBS -N list_ranking
#PBS -l nodes=1:ppn=32:gpus=1
#PBS -j oe
#PBS -q gpu
source /public/software/profile.d/mpi_openmpi-intel-2.1.2.sh
source /public/software/profile.d/cuda10.0.sh
cd $PBS_O_WORKDIR
gcc gen_list_ranking_data.c -o gen_list_ranking_data -std=c99
./gen_list_ranking_data > list_ranking_data.txt
gcc list_ranking.c -o list_ranking -std=c99
time ./list_ranking < list_ranking_data.txt > list_ranking_out.txt
mpicc list_ranking_mpi.c -o list_ranking_mpi -std=c99
time mpiexec -machinefile $PBS_NODEFILE ./list_ranking_mpi < list_ranking_data.txt > list_ranking_mpi_out.txt
nvcc list_ranking_cuda.cu -o list_ranking_cuda
time ./list_ranking_cuda < list_ranking_data.txt > list_ranking_cuda_out.txt
运行时间list_ranking.o17502
自上而下分别是串行算法的时间（1.977s），MPI 算法的时间（5.976s），CUDA 算法的时间（3.246s）。可以看到，这里实现的两个并行算法都在不过分增加并行开销的情况下增加了串行算法的可扩展性。
这个问题中并行算法是不可能快于同等情况下的串行算法的，因为串行算法的时间复杂度是$O(n)$，而读入这个向量的时间复杂度已经是$O(n)$了，并行算法无论如何都需要先读入向量，时间复杂度上不可能优于$O(n)$。
real 0m1.977s
user 0m1.710s
sys 0m0.176s
real 0m5.976s
user 1m20.123s
sys 0m32.837s
real 0m3.246s
user 0m2.251s
sys 0m0.806s
构造生成数据gen_list_ranking_data.c
此处我使用了老师提供的数据生成器。
#include <stdio.h>
#include <stdlib.h>
int *gen_linked_list_1(unsigned int N)
{
int *list = NULL;
if (NULL != list)
{
free(list);
list = NULL;
}
if (0 == N)
{
printf("N is 0, exit\n");
exit(-1);
}
list = (int *)malloc(N * sizeof(int));
if (NULL == list)
{
printf("Can not allocate memory for output array\n");
exit(-1);
}
int i;
for (i = 0; i < N; i++)
list[i] = i - 1;
return list;
}
//i和j的后继元素交换位置
void swap(int *list, int i, int j)
{
if (i < 0 || j < 0 || i == j)
return;
int p = list[i]; //保存i后继元素下标p
int q = list[j]; //保存j后继元素下标q
if (p == -1 || q == -1)
return; //如果有一个没有后继元素
int pnext = list[p]; //保存p的后继元素下标
int qnext = list[q]; //保存q的后继元素下标
//i,j的后继元素交换位置
if (p == j)
{ //j是i的后继
list[i] = q;
list[j] = list[q];
list[q] = j;
}
else if (i == q)
{ //i是j的后继
list[j] = p;
list[i] = list[p];
list[p] = i;
}
else
{
list[i] = q; //i的后继改为q
list[j] = p; //j的后继改为p
list[p] = qnext; //p的后继元素改为原来q的后继
list[q] = pnext; //q的后继元素改为原来p的后继
}
}
int *gen_linked_list_2(unsigned int N)
{
int *list;
list = gen_linked_list_1(N);
int p = N / 5;
int i, temp, k;
for (i = 0; i < N; i += 2)
{
int k = (i + i + p) % N;
swap(list, i, k);
}
return list;
}
int main()
{
int N = 10000000;
int *qq = NULL;
int i;
// qq=gen_linked_list_1(N);
// printf("\nhere is the list\n");
// for(i=0; i<N; i++)
// printf("%3d ", qq[i]);
// printf("\n");
// free(qq);
qq = gen_linked_list_2(N);
//printf("\nhere is the new list\n");
printf("%3d ", N); //输出链表元素个数
printf("\n");
for (i = 0; i < N; i++) //输出链表全部元素
printf("%3d ", qq[i]);
printf("\n");
free(qq);
return 0;
}
终端执行下述指令，得到用于测试的数据list_ranking_data.txt。
gcc gen_list_ranking_data.c -o gen_list_ranking_data -std=c99
./gen_list_ranking_data > list_ranking_data.txt
list_ranking_data.txt有着如下的结构，这里只放出向量的前四个元素。
10000000
-1 4000014 2000003 4000030 ...
串行扫描算法
最简单的串行算法即扫描整个列表。解释一下整个算法的过程：
- 读入向量
v - 根据
v计算每个元素的前驱pre，同时找出链表尾last - 从链表尾
last开始顺次访问前驱pre，遍历的顺序就是我们要求的链表序。
#include <stdio.h>
#include <stdlib.h>
int main()
{
int n, last;
scanf("%d", &n);
int *v = (int *)malloc(sizeof(int) * n),
*a = (int *)malloc(sizeof(int) * n),
*pre = (int *)malloc(sizeof(int) * n);
for (int i = 0; i < n; ++i)
scanf("%d", &v[i]);
for (int i = 0; i < n; ++i)
{
if (v[i] < 0)
last = i;
else
pre[v[i]] = i;
}
for (int u = last, i = n - 1; i >= 0; --i, u = pre[u])
a[u] = i;
for (int i = 0; i < n; ++i)
printf("%d ", a[i]);
free(pre);
free(a);
free(v);
}
终端执行下属指令，可以计时执行并将结果写入list_ranking_out.txt。
gcc list_ranking.c -o list_ranking -std=c99
time ./list_ranking < list_ranking_data.txt > list_ranking_out.txt
list_ranking_out.txt前四个元素为如下的结构，手动检查第一个元素的 ranking 是9999999，是正确的。
9999999 5999984 9999997 5999968 ...
直接串行扫描算法修改出来的并行算法不会有理想的加速效果。理由是，下面这段代码中，重复执行u = pre[u]从链表尾开始依次向前找前驱节点，这个过程是串行的，有循环依赖，不可以直接划分到多进程并行执行。
for (int u = last, i = n - 1; i >= 0; --i, u = pre[u])
a[u] = i;
MPI 实现的 Wyllie 算法list_ranking_mpi.c
即使用 point jumping（点倍增）的方法来计算各个 rank，是最简单的并行算法。
Wyllie 算法原理：
for i in range(n): # 初始化，其中v代表初始向量，即每个节点后继 a[i] = v[i]>=0 # a代表每个节点 while exists(i => v[i]>=0): # 此部分并行执行，每轮迭代要同步 a[i] += a[v[i]] v[i] = v[v[i]]
如上面的伪代码，初始化过程和每轮迭代都是可以划分到多进程执行的。当然，这里得到的a是到链表尾部的距离，用链表长度相减就能得到我们所需要的序号。
#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
int main(int argc, char **argv)
{
int comSize, comRank, n;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &comSize);
MPI_Comm_rank(MPI_COMM_WORLD, &comRank);
if (!comRank)
scanf("%d", &n);
MPI_Bcast(
&n,
1,
MPI_INT,
0,
MPI_COMM_WORLD);
int myLen = (n + comSize - 1) / comSize,
myBeg = comRank * myLen,
*myV = (int *)malloc(sizeof(int) * myLen),
*myA = (int *)malloc(sizeof(int) * myLen),
*v = (int *)malloc(sizeof(int) * n),
*a = (int *)malloc(sizeof(int) * n);
if (!comRank)
for (int i = 0; i < n; ++i)
scanf("%d", &v[i]);
MPI_Bcast(
v,
n,
MPI_INT,
0,
MPI_COMM_WORLD);
for (int i = 0; i < myLen; ++i)
{
myV[i] = v[i + myBeg];
myA[i] = myV[i] >= 0;
}
MPI_Allgather(
myA,
myLen,
MPI_INT,
a,
myLen,
MPI_INT,
MPI_COMM_WORLD);
for (int j = 1; j < n; j <<= 1)
{
for (int i = 0; i < myLen; ++i)
if (v[i + myBeg] >= 0)
{
myA[i] = a[i + myBeg] + a[v[i + myBeg]];
myV[i] = v[v[i + myBeg]];
}
MPI_Allgather(
myA,
myLen,
MPI_INT,
a,
myLen,
MPI_INT,
MPI_COMM_WORLD);
MPI_Allgather(
myV,
myLen,
MPI_INT,
v,
myLen,
MPI_INT,
MPI_COMM_WORLD);
}
if (!comRank)
for (int i = 0; i < n; ++i)
printf("%d ", n - 1 - a[i]);
free(myV);
free(myA);
free(v);
free(a);
MPI_Finalize();
}
终端执行下述指令，可以计时执行并将结果写入list_ranking_mpi_out.txt。
mpicc list_ranking_mpi.c -o list_ranking_mpi -std=c99
time mpiexec -machinefile $PBS_NODEFILE ./list_ranking_mpi < list_ranking_data.txt > list_ranking_mpi_out.txt
list_ranking_mpi_out.txt前四个元素为如下的结构，手动检查第一个元素的 ranking 是9999999，是正确的；list_ranking_mpi_out.txt与list_ranking_out.txt完全相同。
9999999 5999984 9999997 5999968 ...
CUDA 实现的 Wyllie 算法list_ranking_cuda.cu
另一种并行链表求序的 SMP 算法[Helman and JaJa, 1999]需要提前知道链表中等分点的分布，不是很适合本例，因此用 CUDA 实现的仍然是 Wyllie 算法。
相比于 MPI 实现的版本，这里没有了节点广播 Bcast 和全聚集 Allgather 的开销，因此时间上快了非常多。然而，仍然需要考虑的是，核函数 work 对v[v[idx]]和a[v[idx]]的访问没有对齐，这是算法的瓶颈所在。
#include <stdio.h>
#include <stdlib.h>
void __global__ init(int n, int *v, int *a)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
a[idx] = v[idx] >= 0;
}
void __global__ work(int n, int *v, int *a, int *new_v, int *new_a)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (v[idx] >= 0)
{
new_a[idx] = a[idx] + a[v[idx]];
new_v[idx] = v[v[idx]];
}
}
int main()
{
int n;
scanf("%d", &n);
int *v = (int *)malloc(sizeof(int) * n),
*d_v,
*d_v1,
*d_a,
*d_a1,
block = 32, grid = n / block;
for (int i = 0; i < n; ++i)
scanf("%d", &v[i]);
cudaMalloc(&d_a, sizeof(int) * n);
cudaMalloc(&d_a1, sizeof(int) * n);
cudaMalloc(&d_v, sizeof(int) * n);
cudaMalloc(&d_v1, sizeof(int) * n);
cudaMemcpy(d_v, v, sizeof(int) * n, cudaMemcpyHostToDevice);
init<<<grid, block>>>(n, d_v, d_a);
cudaDeviceSynchronize();
for (int j = 1; j < n; j <<= 1)
{
work<<<grid, block>>>(n, d_v, d_a, d_v1, d_a1);
cudaDeviceSynchronize();
int *tmp;
tmp = d_v, d_v = d_v1, d_v1 = tmp;
tmp = d_a, d_a = d_a1, d_a1 = tmp;
}
cudaMemcpy(v, d_a, sizeof(int) * n, cudaMemcpyDeviceToHost);
for (int i = 0; i < n; ++i)
printf("%d ", n - 1 - v[i]);
cudaFree(d_v);
cudaFree(d_a);
cudaFree(d_v1);
cudaFree(d_a1);
free(v);
}
终端执行下述指令，可以计时执行并将结果写入list_ranking_cuda_out.txt。
nvcc list_ranking_cuda.cu -o list_ranking_cuda
time ./list_ranking_cuda < list_ranking_data.txt > list_ranking_cuda_out.txt
list_ranking_cuda_out.txt前四个元素为如下的结构，手动检查第一个元素的 ranking 是9999999，是正确的；list_ranking_cuda_out.txt与list_ranking_out.txt完全相同。
9999999 5999984 9999997 5999968 ...
参考文献
- http://staff.ustc.edu.cn/~csli/graduate/algorithms/book6/chap30.htm
- https://en.wikipedia.org/wiki/Pointer_jumping
- http://www.cs.cmu.edu/~scandal/alg/listrank.html
- https://www.cs.cmu.edu/~glmiller/Publications/Papers/ReMiMo93.pdf
- http://cdn.iiit.ac.in/cdn/cstar.iiit.ac.in/~kkishore/ics152-rehman.pdf
- D. Bader, G. Cong, and J. Feo. On the Architectural Requirements for Efficient Execution of Graph Algorithms. International Conference on Parallel Processing (ICPP), 2005, pages 547-556, June 2005.
- D. A. Bader, V. Agarwal, and K. Madduri. On the Design and Analysis of Irregular Algorithms on the Cell Processor: A Case Study of List Ranking. In 21st IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 1-10. IEEE, 2007.
- D. A. Bader and G. Cong. A Fast, Parallel Spanning Tree Algorithm for Symmetric Multiprocessors (SMPs). Journal of Parallel and Distributed Computing, 65(9):994 - 1006, 2005.
- https://link.springer.com/content/pdf/10.1007%2FBFb0056600.pdf