双向Dijkstra算法的MPI+Openmp实现

Reference

Readme

运行方法

  • mpic++ -o DiDijkstra_heap DiDijkstra_heap.cpp -fopenmp -mcmodel=large -O3
  • time mpiexec -n 3 ./DiDijkstra_heap 20001, 必须要指定图的点的数目, 最大支持20001, 必须设置mpi进程数目为3, 因为本程序特别为3主机情况设计.

注意事项

  • 在所有主机与执行文件同目录必须有名为input.txt的图数据, 数据形式为. 特别注意, 由于手写了输入挂, 所以文件结尾必须为数字.

    1
    2
    3
    4
    0 9 2
    0 10 5
    0 14 9
    ...
  • 多次运行会越来越快, 可以试一下(因为mmap)

  • 默认读入为有向图

Intro

本次实验环境为3台核心数为4的服务器, 要求使用openmp+mpi进行编程. 问题规模为2w7000w边. 在进行简单测试后, 我认为在这个规模下问题主要的耗时为:

  • 文件IO耗时
  • MPI通信耗时
  • 如果MPI操作为同步操作, 会非常耗时, 所以需要改进算法的通信为异步通信.

我使用双向Dijkstra方法, 这个经过我特别设计的双向算法的特点是保证最优解, 好处是算法运行时不需要同步通信, 而且搜索范围直观上非常小, 速度极快. 而且由于这个异步特性, 可以自动均衡两台电脑上的运行时间, 保证两台worker电脑从开始到结束都是在充分运行的. 当然代价是master基本在空转, 不过我做这个项目只是为了实现我的想法, 现在看来效果似乎还不错, 不过, 由于问题规模太小, 算法威力难以体现, 如果这个图非常大, 边非常多, 我的算法在并行实现上就会十分有优势, 因为它直接减少了搜索的点的数目(在理想情况下).

双向Dijkstra算法

算法流程

  1. 分别从起点和终点开始运行dijkstra算法
  2. 当两个Dijkstra算法找到了一个重叠点, 记为(*)
  3. 此时两个Dijkstra算法接着运行, 但不再利用边更新距离数组.
  4. 如果找到新的重叠点, 则需要与(*)比较, 选出最优的那个作为新的(*)
  5. 最终产生的*必定为最优解

这一点只能通过修改算法实现, 可以使用双向dj(因为发送数据可以不要同步), 一个从头开始, 一个从尾开始, 每一次都将自己选出的点发给master, 当选出的点重复时, 就表明找到一条可能是最短路的路径. 这时候master告诉worker, worker此时不再更新dist, 而是在已经更新过的点中不断产生nearest_node, 并看看有没有更好的重叠点, 等两个dj都跑完, 这时候就找到了最优解.

算法正确性

分析, 当找到第一个重叠点时, 现在有的点分三类

  • 一类点是已经被选择为nearest_node的点.
  • 二类点是已经被更新过dist选为过的, 这一类的点在之后作为算法选出的新的nearest_node的过程中, 到起点的距离是不断增加的.
  • 三类点是尚未被访问过的, 连dist都是正无穷.

那么, 现在新的可能的最短路只可能在新选出的二类点与已有的一类点重叠产生. 我们接下来将证明这个结论. 因为如果是这样, 我们只需要遍历完此时剩下的二类点, 看看它们有没有产生一类点与二类点的重叠. 如果有, 再判断它们是否构成更好的解.

  • (1,1): 一类点和一类点的重叠点不可能为更好的解, 因为我们是在找到第一个一类点与一类点重叠时考虑的.
  • (1,2): 这个是我们算法想要处理的剩余情况, 所以不用考虑.
  • (2,2), (2,3), (3,3): 这些情况下, 点分别到起点和终点的最短距离都大于(*)到起点和终点的最短距离, 所以也不会是更好的解
  • (1,3): 那么, 现在的问题是, 有没有可能, 一类点与三类点重叠形成的路径比最开始找到的路径更短?答案是没有, 因为, 如果有这么一个重叠点, 它所在的路径, 必定要经过(2,2)或者(2,3)(如果是经过(1,2)的情况, 则必定可以被我们的算法搜索到, 所以不用考虑), 而如果再假设(1,3)重叠点是最优解, 由最优解路径上的点到起点和终点的距离都是最短距离, 那么考虑这样一条路径上的(2,2)点, 则它分别到起点和终点的距离又是最短路径, 由(2,2)点分别到起点和终点的最短路径都要大于(*)分别到起点和终点的最短路径. 所以(1,3)重叠点不可能比(*)更优.

所以, 我们只需要让worker在知道自己选出的点重叠后不再更新dist(不增加新的二类点)就可以了.

至于(1,3)的最短路情况为什么必定要经过(1,2)或(2,2)或(2,3), 这个可以通过画图轻松看出.

以上分析是很粗糙的, 可能结论有错误.

优化手段

异步通信代替同步

因为我设计的算法的特殊性, 可以在发送当前选出的点的时候不需要同步

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#ifdef TIMER_MPI
timespec time1, time2;
clock_gettime(CLOCK_REALTIME, &time1);
#endif
// MPI send curnode
{
worker_tuples[cnt][0] = nearest_vex;
if (nearest_vex != -1) worker_tuples[cnt][1] = worker_dist[nearest_vex];
MPI_Isend(worker_tuples[cnt], 2, MPI_INT, 0, CurNode, comm, &request);
++cnt;
}
#ifdef TIMER_MPI
clock_gettime(CLOCK_REALTIME, &time2);
timespec diff(time_diff(time1, time2));
c_cnt += diff;
#endif

worker间的异步->自动负载均衡

如果进行测试, 会发现如果一个workerIO用时很久, 那么它Dijkstra用时就会比较小, 这一点可以反复测试一下:

1
2
3
4
5
6
time mpiexec -n 3 ./DiDijkstra_heap 20001
reader 1 use: 4.473134170 sec
reader 2 use: 4.711733391 sec
worker 1 use: 0.241378561 sec
worker 2 use: 0.2210793 sec
0 17665 20000

这是因为我在两个workerIO后没有做同步. 而是让它们里面跑Dijkstra, 因为双向Dijkstra不需要同步. 而在双向Dijkstra中算法运行的临界是第一次有重叠点, 那么, 由于先IO完的先搜索, 那它搜索的范围就会比另一个worker搜索的范围大, 就自动减轻了另一个worker的工作量. (然而这个数据集太感人了, 最短路居然是两跳, 很难体现出这个优越性)

不使用memset和malloc

虽然malloc开销很小, 但是malloc的空间需要memset, 而memset的开销很大. 所以采用静态全局数组代替动态内存分配, 这样有个好处, 数组元素默认为0, 不需要memset.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
int n;
int table[N][N][2];
int table_cnt[N];

int master_cnt[N];
int master_dist[N];
int master_res1[N];
int master_res2[N];
void master(int n, int comm);

bool worker_visited[N];
int worker_from[N];
int worker_dist[N];
int worker_tuples[N][2];
Heap<Node, 100*N> worker_heap;
void worker(int my_rank, int n, int comm, int v1, int num_threads);

邻接表代替邻接矩阵

这样做在IO时增加的工作量不大, 但却能很好地减少Dijkstra遍历边的常数, 特别是在这个情景下边相对稀疏.

1
2
3
index = table_cnt[from]++; \\ 记录邻接表这一项有多少已用过的位置
table[from][index][0] = to;
table[from][index][1] = weight;

减少了这里的常数

1
2
3
4
5
6
7
8
9
10
11
12
if (flag != -1) {
#pragma omp parallel for
for (int i = 0; i < table_cnt[nearest_vex]; i++) {
int to = table[nearest_vex][i][0];
int edge_weight = table[nearest_vex][i][1];
int dist = worker_dist[nearest_vex] + edge_weight;
if (worker_dist[to] > dist) {
worker_dist[to] = dist;
worker_heap.push(Node(nearest_vex, to, dist));
}
}
}

堆优化

在找最近点的时候, 可以使用堆优化

1
2
3
4
5
6
7
8
9
nearest_vex = -1; min_dist = 0x7f7f7f7f;
while (!worker_heap.is_empty()) {
Node node = worker_heap.top();
worker_heap.pop();
if (worker_visited[node.to]) continue;
worker_from[node.to] = node.from;
nearest_vex = node.to;
break;
}

这个堆是我自己手写的.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
template<typename T, size_t Size>
class Heap
{
public:
Heap()
: mSize(0)
{

}

void push(const T& data)
{
int i = mSize++; mTree[mSize] = data;
while (i != 0 && data > mTree[(i - 1) / 2]) {
mTree[i] = mTree[(i - 1) / 2];
i = (i - 1) / 2;
}
mTree[i] = data;

return;
}

void pop()
{
if (mSize == 0) return;

T back_elem = mTree[--mSize];

int i = 0; int ci = i * 2 + 1;
while (ci < mSize - 1) {
if (mTree[ci] < mTree[ci + 1]) ci++;

if (mTree[ci] > back_elem) {
mTree[i] = mTree[ci];
i = ci;
ci = ci * 2 + 1;
}
else
break;
}

mTree[i] = back_elem;

return;
}

T top() const
{
return mTree[0];
}

bool is_empty() const
{
return mSize == 0;
}

private:
size_t mSize;
T mTree[Size];
};

文件IO改进

在我测试的情况来看, IO和处理用时90%. 其中特别是缺页读硬盘的部分耗时极高.

在linux中大文件快速读取可以使用freadmmap, 经测试, mmap在此问题下较快. mmap将文件映射到内存, 利用缺页异常读取文件, 可以减少内存拷贝. 这个函数在我的电脑上, 不开O3单进程运行费时20s.

下面比较几种IO手段, 并选择最优. 最终选择方案是单线程的mmap加手写输入挂.

fscanf

1
2
3
4
5
6
7
8
9
void file_read_fscanf()
{
FILE* fd = fopen(file_name, "r");

int from, to, weight;
while (fscanf(fd, "%d %d %d\n", &from, &to, &weight) != EOF) {
table[from * n + to] = weight;
}
}
1
用时 21.897 s

mmap

1
2
单线程: 20.518秒
多线程: 38秒

https://www.byvoid.com/zhs/blog/fast-readfile

值得注意的是, 程序运行过一次后, 由于文件已经被加载入内存, 再次读不会发送缺页, 所以会很快, 相当于直接做内存访问, 这时候最快速度是5s. 如果是热读, 多线程更快(1.4s), 但如果是冷读, 单线程更快(20s).

1
2
3
4
5
6
7
8
9
10
11
void file_read_mmap()
{
int fd = open(file_name, O_RDONLY);
int size = lseek(fd, 0, SEEK_END);
char *mbuf = (char *) mmap( NULL, size ,PROT_READ, MAP_SHARED, fd, 0 );
char* end = mbuf + size;

read_worker(mbuf, end);

munmap(mbuf, size);
}

我们来看看read_worker, 这个函数是特别针对这个问题的输入形式编写的, 一次读取一行. 利用宏展开自动内联, 利用位运算和前置自加运算符加快速度.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#define GET_INT(res, ch, mbuf) \
do{\
res = 0;\
ch = *mbuf;\
while ('0'<=ch && ch<='9') {\
res = (res<<3) + (res<<1) + (ch^'0');\
ch = *(++mbuf);\
}\
++mbuf;\
}while(false);\

void read_worker(char* mbuf, char* end)
{
int from, to, weight;
char ch;
while (mbuf < end) {
GET_INT(from, ch, mbuf);
GET_INT(to, ch, mbuf);
GET_INT(weight, ch, mbuf);

table[from * n + to] = weight;
}
}
文章目录
  1. 1. Reference
  2. 2. Readme
  3. 3. Intro
  4. 4. 双向Dijkstra算法
    1. 4.1. 算法流程
    2. 4.2. 算法正确性
  5. 5. 优化手段
    1. 5.1. 异步通信代替同步
    2. 5.2. worker间的异步->自动负载均衡
    3. 5.3. 不使用memset和malloc
    4. 5.4. 邻接表代替邻接矩阵
    5. 5.5. 堆优化
    6. 5.6. 文件IO改进
      1. 5.6.1. fscanf
      2. 5.6.2. mmap
|