OpenMP 旅行商问题，静态调度
- 前端
- Stella981
- 0
▶ 《并行程序设计导论》第六章中讨论了旅行商，分别使用了 MPI，Pthreads，OpenMP 来进行实现，这里是 OpenMP 的代码，分为静态调度（每个线程分分配等量的搜索人物）和动态调度（每个线程分配不等量的任务，每当有线程完成自己的任务后，向其他线程请求新的子任务）
● 静态调度代码
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <omp.h>
5
6 #define DEBUG
7 #define City_count(tour) (tour->count)
8 #define Tour_cost(tour) (tour->cost)
9 #define Last_city(tour) (tour->cityList[(tour->count)-1]) // tour 中最后一个城市的编号
10 #define Tour_city(tour,i) (tour->cityList[(i)]) // tour 中第 i 个城市的编号
11 #define Cost(city1, city2) (digraph[city1 * n + city2]) // tour 中两个城市之间的代价
12 #define Queue_elt(queue,i) (queue->list[(queue->head + (i)) % queue->spaceAlloc])
13
14 typedef int city_t; // 城市编号
15 typedef int cost_t; // 代价，可以为浮点数
16
17 typedef struct
18 {
19 int count; // 已经走过的城市数
20 city_t* cityList; // 已经走过的城市列表
21 cost_t cost; // 当前代价
22 } tour_struct;
23
24 typedef struct
25 {
26 int spaceAlloc; // 栈占据的空间大小（最大大小）
27 int stackSize; // 栈的实际大小
28 tour_struct** list;
29 } stack_struct;
30
31 typedef struct
32 {
33 int spaceAlloc; // 队列占据的空间大小（最大大小）
34 int head; // 队头
35 int tail; // 队尾
36 bool full; // 是否队满
37 tour_struct** list;
38 } queue_struct;
39
40 const int INFINITY = 1000000; // 最大代价，用于初始化 tour
41 const int NO_CITY = -1; // 标记 tour->cityList 中尚未使用的部分
42 const int MAX_STRING = 1000; // 输出字符串最大长度
43 int n; // 城市总数
44 int nThread; // 线程数
45 city_t home_town = 0; // 起点和终点的编号
46 cost_t *digraph; // 读取的邻接表数据
47 tour_struct *best_tour; // 保存最优轨迹
48 queue_struct *queue; // 搜索队列
49 int queueSize; // 搜索队的大小
50 int initialTourCount; // 首次分配给各线程的搜索队列数量
51 tour_struct* allocTour(stack_struct* avail);// 必需的声明，不然该函数与栈操作函数混在一起
52
53 // 基本函数
54 void usage(char* prog_name)// 提示信息
55 {
56 fprintf(stderr, "\n<Usage> %s <nThread> <digraph file>\n", prog_name);
57 exit(0);
58 }
59
60 void readDigraph(FILE* digraph_file)// 从文件读取邻接表
61 {
62 int i, j;
63 fscanf_s(digraph_file, "%d", &n);// 第一行数字时城市数量
64 if (n <= 0)
65 {
66 fprintf(stderr, "\n<readDigraph> Number of city = %d\n", n);
67 exit(1);
68 }
69 digraph = (cost_t*)malloc(n * n * sizeof(cost_t));
70 for (i = 0; i < n; i++)
71 {
72 for (j = 0; j < n; j++)
73 {
74 fscanf_s(digraph_file, "%d", &digraph[i * n + j]);
75 if (i == j && digraph[i * n + j] != 0)
76 {
77 fprintf(stderr, "\n<readDigraph> Diagonal element [%d, %d] = %d\n", i, j, digraph[i * n + j]);
78 exit(2);
79 }
80 else if (i != j && digraph[i * n + j] <= 0)
81 {
82 fprintf(stderr, "\n<readDigraph> Off-diagonal element [%d, %d] = %d\n", i, j, digraph[i * n + j]);
83 exit(3);
84 }
85 }
86 }
87 #ifdef DEBUG
88 printf("\n\t<readDigraph> Read map\n");
89 #endif
90 }
91
92 void printDigraph(void)// 打印邻接表
93 {
94 int i, j;
95 printf("\n\t<printDigraph> Order = %d\nMatrix = \n", n);
96 for (i = 0; i < n; i++)
97 {
98 for (j = 0; j < n; j++)
99 printf("%2d ", digraph[i * n + j]);
100 printf("\n");
101 }
102 printf("\n");
103 }
104
105 void initeTour(tour_struct* tour, cost_t cost)// 初始化搜索轨迹，要求给出代价
106 {
107 int i;
108 tour->cityList[0] = 0;
109 for (i = 1; i <= n; tour->cityList[i++] = NO_CITY);
110 tour->cost = cost;
111 tour->count = 1;
112 }
113
114 void copyTour(tour_struct* tour1, tour_struct* tour2)// 搜索轨迹 tour1 拷贝给 tour2
115 {
116 memcpy(tour2->cityList, tour1->cityList, (n + 1) * sizeof(city_t));
117 tour2->count = tour1->count;
118 tour2->cost = tour1->cost;
119 }
120
121 void printTour(int my_rank, tour_struct* tour, char* title)// 输出搜索轨迹，注意有两种输出
122 {
123 int i;
124 char string[MAX_STRING];
125 if (my_rank >= 0)// 各线程汇报
126 sprintf_s(string, "\n\t<printTour> Rank%2d, %s %p: ", my_rank, title, tour);
127 else // 输出当前最佳回路
128 sprintf_s(string, "\n\t<printTour> %s: ", title);
129 for (i = 0; i < City_count(tour); i++)
130 sprintf_s(string + strlen(string), 10, "%d ", Tour_city(tour, i));
131 printf("%s\n", string);
132 }
133
134 // 栈相关
135 stack_struct* initeStack(void)// 初始化栈
136 {
137 int i;
138 stack_struct* stack = (stack_struct*)malloc(sizeof(stack_struct));
139 stack->list = (tour_struct**)malloc(n * n * sizeof(tour_struct*));
140 for (i = 0; i < n * n; stack->list[i++] = NULL);
141 stack->stackSize = 0;
142 stack->spaceAlloc = n * n;
143 return stack;
144 }
145
146 void pushStack(stack_struct* stack, tour_struct* tour)// 非拷贝压栈
147 {
148 if (stack->stackSize == stack->spaceAlloc)
149 {
150 fprintf(stderr, "\n<pushStack> Overflow\n");
151 free(tour->cityList);
152 free(tour);
153 }
154 else
155 {
156 # ifdef DEBUG
157 printf("\n\t<pushStack> StackSize = %d, tour at %p, tour->cityList at %p\n", stack->stackSize, tour, tour->cityList);
158 printTour(-1, tour, "pushStack");
159 printf("\n");
160 # endif
161 stack->list[stack->stackSize] = tour;
162 (stack->stackSize)++;
163 }
164 }
165
166 void pushCopyStack(stack_struct* stack, tour_struct* tour, stack_struct* avail)// 拷贝压栈
167 {
168 tour_struct* tmp;
169 if (stack->stackSize == stack->spaceAlloc)
170 {
171 fprintf(stderr, "\n<pushCopyStack> Overflow\n");
172 exit(-1);
173 }
174 tmp = allocTour(avail);
175 copyTour(tour, tmp);
176 stack->list[stack->stackSize] = tmp;
177 (stack->stackSize)++;
178 }
179
180 tour_struct* popStack(stack_struct* stack)// 出栈
181 {
182 tour_struct* tmp;
183 if (stack->stackSize == 0)
184 {
185 fprintf(stderr, "\n<popStack> Empty\n");
186 exit(-1);
187 }
188 tmp = stack->list[stack->stackSize - 1];
189 stack->list[stack->stackSize - 1] = NULL;
190 (stack->stackSize)--;
191 return tmp;
192 }
193
194 int isEmptyStack(stack_struct* stack)// 判断是否空栈
195 {
196 return stack->stackSize == 0;
197 }
198
199 void freeStack(stack_struct* stack)// 清空栈
200 {
201 int i;
202 for (i = 0; i < stack->stackSize; free(stack->list[i]->cityList), free(stack->list[i]), i++);
203 free(stack->list);
204 free(stack);
205 }
206
207 void printStack(stack_struct* stack, int my_rank, char title[])// 打印栈
208 {
209 char string[MAX_STRING];
210 int i, j;
211 printf("\n\t<printStack> Rank%2d, %s\n", my_rank, title);
212 for (i = 0; i < stack->stackSize; i++)
213 {
214 sprintf_s(string, 10, "%d> ", i);
215 for (j = 0; j < stack->list[i]->count; j++)
216 sprintf_s(string + strlen(string), 10, "%d ", stack->list[i]->cityList[j]);
217 printf("%s\n", string);
218 }
219 }
220
221 // 队列相关
222 queue_struct* initeQueue(int size)// 初始化队列
223 {
224 queue_struct* new_queue = (queue_struct*)malloc(sizeof(queue_struct));
225 new_queue->list = (tour_struct**)malloc(size * sizeof(tour_struct*));
226 new_queue->spaceAlloc = size;
227 new_queue->head = new_queue->tail = new_queue->full = false;
228 return new_queue;
229 }
230
231 int isEmptyQueue(queue_struct* queue)// 判断是否队空
232 {
233 return !queue->full && queue->head == queue->tail;// 空队要求 queue->full == false 且头尾指针相等
234 }
235
236 tour_struct* deQueue(queue_struct* queue)// 出队
237 {
238 tour_struct* tmp;
239 if (isEmptyQueue(queue))
240 {
241 fprintf(stderr, "\n<deQueue> Empty queue\n");
242 exit(-1);
243 }
244 tmp = queue->list[queue->head];
245 queue->head = (queue->head + 1) % queue->spaceAlloc;
246 return tmp;
247 }
248
249 void enQueue(queue_struct* queue, tour_struct* tour)// 入队
250 {
251 tour_struct* tmp;
252 if (queue->full == true)
253 {
254 fprintf(stderr, "\n<enQueue> Overflow\n");
255 exit(-1);
256 }
257 tmp = allocTour(NULL);
258 copyTour(tour, tmp);
259 queue->list[queue->tail] = tmp;
260 queue->tail = (queue->tail + 1) % queue->spaceAlloc;
261 if (queue->tail == queue->head)
262 queue->full = true;
263 }
264
265 void freeQueue(queue_struct* queue)// 清空队列
266 {
267 free(queue->list);
268 free(queue);
269 }
270
271 void printQueue(queue_struct* queue, int my_rank, char title[])// 打印队列
272 {
273 char string[MAX_STRING];
274 int i, j;
275 printf("\n\t<printQueue> Rank%2d > %s\n", my_rank, title);
276 for (i = queue->head; i != queue->tail; i = (i + 1) % queue->spaceAlloc)
277 {
278 sprintf_s(string, "%d> %p = ", i, queue->list[i]);
279 for (j = 0; j < queue->list[i]->count; j++)
280 sprintf_s(string + strlen(string), 10, "%d ", queue->list[i]->cityList[j]);
281 printf("%s\n", string);
282 }
283 }
284
285 // Tour 相关
286 tour_struct* allocTour(stack_struct* avail)// 生成一个搜索轨迹
287 {
288 tour_struct* tmp;
289 if (avail == NULL || isEmptyStack(avail))
290 {
291 tmp = (tour_struct*)malloc(sizeof(tour_struct));
292 tmp->cityList = (city_t*)malloc((n + 1) * sizeof(city_t));
293 return tmp;
294 }
295 return popStack(avail);
296 }
297
298 int visited(tour_struct* tour, city_t city)// 判断一个轨迹中是否经过了城市 city
299 {
300 for (int i = 0; i < City_count(tour); i++)
301 {
302 if (Tour_city(tour, i) == city)
303 return true;
304 }
305 return false;
306 }
307
308 void addCity(tour_struct* tour, city_t new_city)// 将一个城市添加到 tour
309 {
310 city_t old_last_city = Last_city(tour); // 从原来的最后一个城市接出一条路径
311 tour->cityList[tour->count] = new_city;
312 (tour->count)++;
313 tour->cost += Cost(old_last_city, new_city);
314 }
315
316 void removeLastCity(tour_struct* tour)// 去掉 tour 中最后一个城市
317 {
318 city_t old_last_city = Last_city(tour), new_last_city;
319 tour->cityList[tour->count - 1] = NO_CITY;
320 (tour->count)--;
321 new_last_city = Last_city(tour);
322 tour->cost -= Cost(new_last_city, old_last_city);
323 }
324
325 int isBestTour(tour_struct* tour)// 判断当前轨迹是否最佳
326 {
327 cost_t cost_so_far = Tour_cost(tour);
328 city_t last_city = Last_city(tour);
329 return cost_so_far + Cost(last_city, home_town) < Tour_cost(best_tour);// 当前轨迹的代价加上最后一个城市回家的代价是否更优
330 }
331
332 int isFeasible(tour_struct* tour, city_t city)// 判断 tour 的最后一个城市是否可以前往城市 city
333 {
334 return !visited(tour, city) && Tour_cost(tour) + Cost(Last_city(tour), city) < Tour_cost(best_tour);// 要求没去过该城市且代价不超出当前最优回路代价
335 }
336
337 void updateBestTour(tour_struct* tour)
338 {
339 if (isBestTour(tour))// 再次检查是否最优，防止两次检查之间其他线程更新了最优轨迹
340 {
341 copyTour(tour, best_tour);
342 addCity(best_tour, home_town);
343 }
344 }
345
346 void freeTour(tour_struct* tour, stack_struct* avail)
347 {
348 if (avail == NULL) // 普通释放
349 {
350 free(tour->cityList);
351 free(tour);
352 }
353 else // 将结点压回栈中
354 pushStack(avail, tour);
355 }
356
357 void distributeTour(int my_rank, int* my_first_tour_p, int* my_last_tour_p)// 分配初始搜索任务，使用块划分
358 {
359 const int quotient = initialTourCount / nThread, remainder = initialTourCount % nThread;// 平均每线程 quotient 个，余数部分由靠前的线程分担
360 int my_count;
361 if (my_rank < remainder)// 靠前的线程，分担余数
362 {
363 my_count = quotient + 1;
364 *my_first_tour_p = my_rank*my_count;
365 }
366 else // 靠后的线程，直接分配，注意偏移
367 {
368 my_count = quotient;
369 *my_first_tour_p = my_rank*my_count + remainder;
370 }
371 *my_last_tour_p = *my_first_tour_p + my_count - 1;
372 }
373
374 // 核心计算函数附属
375 inline long long factorial(int k)// 计算 k 的阶乘
376 {
377 long long tmp = 1;
378 int i;
379 for (i = 2; i <= k; tmp *= i++);
380 return tmp;
381 }
382
383 inline int supQueueSize(void)// 估计搜索需要的队列长度上限，并判断使用的线程数是否太多，没看懂
384 {
385 int fact, size;
386 for (fact = size = n - 1; size < nThread; size *= ++fact);
387 if (size > factorial(n - 1))
388 {
389 fprintf(stderr, "\n<supQueueSize> Too many thread\n");
390 size = 0;
391 }
392 return size;
393 }
394
395 void buildInitialQueue(void)// 生成初始搜索队列
396 {
397 int currentQueueSize; // 当前队列的长度
398 city_t nbr;
399 tour_struct* tour = allocTour(NULL);
400
401 initeTour(tour, 0); // 仅包含起点，代价为 0
402 queue = initeQueue(2 * queueSize); // 搜索队
403
404 enQueue(queue, tour); // 将起点拷贝入队
405 freeTour(tour, NULL);
406 for (currentQueueSize = 1; queueSize < nThread;)// 向搜索队中添加元素，直到不小于线程数为止，以便分配各线程工作
407 {
408 tour = deQueue(queue); // 从队列中取出一个城市来
409 currentQueueSize--;
410 for (nbr = 1; nbr < n; nbr++)
411 {
412 if (!visited(tour, nbr)) // 取出的城市是新城市
413 {
414 addCity(tour, nbr); // 将该城市添加进 tour
415 enQueue(queue, tour); // 将当前 tour 拷贝加入搜索队列中
416 currentQueueSize++;
417 removeLastCity(tour); // 加入搜索队之后立即从 tour 中清楚刚加入的城市，tour 仅用于搜索队操作
418 }
419 }
420 freeTour(tour, NULL); // 本轮入队结束
421 }
422 initialTourCount = currentQueueSize;// 记录当前搜索队中 tour 的数量
423 # ifdef DEBUG
424 printQueue(queue, 0, "\n\t<buildInitialQueue> Queue Initialized\n");
425 # endif
426 }
427
428 void treePartition(int my_rank, stack_struct* stack)// 分树
429 {
430 int my_first_tour, my_last_tour, i;
431 # pragma omp single // 估计搜索队列长度上限
432 queueSize = supQueueSize();
433 # ifdef DEBUG
434 printf("\n\t<treePartition> Rank%2d> supQueueSize = %d\n", my_rank, queueSize);
435 # endif
436 if (queueSize == 0) // 线程数太多
437 exit(0);
438
439 # pragma omp master // 主线程创建初始搜索队列
440 buildInitialQueue();
441 # pragma omp barrier // 从线程要显式的等待主线程
442
443 distributeTour(my_rank, &my_first_tour, &my_last_tour);// 分配初始搜索任务
444 # ifdef DEBUG
445 printf("\n\t<treePartition> Rank%2d> initialTourCount = %d, first = %d, last = %d\n", my_rank, initialTourCount, my_first_tour, my_last_tour);
446 # endif
447 for (i = my_last_tour; i >= my_first_tour; i--)// 从分配到的任务列中中最后一个向前逐个压入栈中
448 {
449 # ifdef DEBUG
450 printTour(my_rank, Queue_elt(queue, i), "pushStack");
451 # endif
452 pushStack(stack, Queue_elt(queue, i));
453 }
454 # ifdef DEBUG
455 printStack(stack, my_rank, "Task set up");
456 # endif
457
458 }
459
460 // 核心计算函数
461 void treeSearch(void)
462 {
463 int my_rank = omp_get_thread_num();
464 city_t nbr;
465 stack_struct *stack = initeStack(), *avail = initeStack();// 搜索栈和未搜索的栈
466 tour_struct *curr_tour;
467
468 for (treePartition(my_rank, stack); !isEmptyStack(stack);)// 初次分树，然后全部搜索完以前不断循环
469 {
470 curr_tour = popStack(stack);
471 # ifdef DEBUG
472 printTour(my_rank, curr_tour, "popStack");
473 # endif
474 if (City_count(curr_tour) == n) // 当前回路已经搜索了 n 个城市
475 {
476 if (isBestTour(curr_tour)) // 判断是否是最佳回路
477 {
478 # ifdef DEBUG
479 printTour(my_rank, curr_tour, "Best tour");
480 # endif
481 # pragma omp critical
482 updateBestTour(curr_tour);// 更新当前最佳回路
483 }
484 }
485 else // 还没有搜索 n 个城市
486 {
487 for (nbr = n - 1; nbr >= 1; nbr--)
488 {
489 if (isFeasible(curr_tour, nbr))
490 {
491 addCity(curr_tour, nbr);
492 pushCopyStack(stack, curr_tour, avail);
493 removeLastCity(curr_tour);
494 }
495 }
496 }
497 freeTour(curr_tour, avail);
498 }
499 freeStack(stack); // 释放线程栈资源
500 freeStack(avail);
501 # pragma omp barrier // 等待所有线程都完成
502 # pragma omp master
503 freeQueue(queue); // 主线程释放队列资源
504 }
505
506 int main(int argc, char* argv[])
507 {
508 FILE* digraph_file; // 输入文件名
509 double start, finish; // 计时器
510 # ifdef DEBUG
511 argc = 3;
512 argv[1] = "1";
513 argv[2] = "D:\\Code\\并行程序设计导论 - 代码\\ch6\\TSP\\mat_17e";
514 # endif
515
516 if (argc != 3) // 检查输入
517 usage(argv[0]);
518 if ((nThread = strtol(argv[1], NULL, 10)) <= 0)
519 {
520 fprintf(stderr, "\n<main> Error thread number\n");
521 usage(argv[0]);
522 }
523 fopen_s(&digraph_file, argv[2], "r");
524 if (digraph_file == NULL)
525 {
526 fprintf(stderr, "\n<main> Error opening input file\n");
527 usage(argv[0]);
528 }
529
530 readDigraph(digraph_file);// 读取图
531 fclose(digraph_file);
532 # ifdef DEBUG
533 printDigraph();
534 # endif
535
536 best_tour = allocTour(NULL);
537 initeTour(best_tour, INFINITY);
538 # ifdef DEBUG
539 printTour(-1, best_tour, "Best tour");
540 printf("\n\t<main> City count = %d\nCost = %d\n\n", City_count(best_tour), Tour_cost(best_tour));
541 # endif
542
543 start = omp_get_wtime();// 开始计算
544 # pragma omp parallel num_threads(nThread) default(none)
545 treeSearch();
546 finish = omp_get_wtime();
547
548 printTour(-1, best_tour, "Best tour");// 汇报结果
549 printf("\n\t<main> Cost = %d\nElapsed time = %e s\n", best_tour->cost, finish - start);
550
551 free(best_tour->cityList);// 释放资源
552 free(best_tour);
553 free(digraph);
554 return 0;
555 }
● 输入的邻接表图
17
0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1
2 0 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 0 2 1 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 0 2 1 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 0 2 1 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 0 2 1 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 0 2 1 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 0 2 1 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 0 2 1 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 0 2 1 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 0 2 1 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 0 2 1 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 0 2 1 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 1 2
1 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2
2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 0 2
2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0
● 输出结果，因为是指数级的枚举，没有等待程序运行结束
免责申明：本站发布的内容（图片、视频和文字）以转载和分享为主，文章观点不代表本站立场，如涉及侵权请联系站长邮箱：xbc-online@qq.com进行反馈，一经查实，将立刻删除涉嫌侵权内容。
