▶ 《并行程序设计导论》第六章中讨论了旅行商，分别使用了 MPI，Pthreads，OpenMP 来进行实现，这里是 OpenMP 的代码，分为静态调度（每个线程分分配等量的搜索人物）和动态调度（每个线程分配不等量的任务，每当有线程完成自己的任务后，向其他线程请求新的子任务）
● 静态调度代码
1 1 #include <stdio.h> 2 2 #include <stdlib.h> 3 3 #include <string.h> 4 4 #include <omp.h> 5 5 6 6 #define DEBUG 7 7 #define City_count(tour) (tour->count) 8 8 #define Tour_cost(tour) (tour->cost) 9 9 #define Last_city(tour) (tour->cityList[(tour->count)-1]) // tour 中最后一个城市的编号 10 10 #define Tour_city(tour,i) (tour->cityList[(i)]) // tour 中第 i 个城市的编号 11 11 #define Cost(city1, city2) (digraph[city1 * n + city2]) // tour 中两个城市之间的代价 12 12 #define Queue_elt(queue,i) (queue->list[(queue->head + (i)) % queue->spaceAlloc]) 13 13 14 14 typedef int city_t; // 城市编号 15 15 typedef int cost_t; // 代价，可以为浮点数 16 16 17 17 typedef struct 18 18 { 19 19 int count; // 已经走过的城市数 20 20 city_t* cityList; // 已经走过的城市列表 21 21 cost_t cost; // 当前代价 22 22 } tour_struct; 23 23 24 24 typedef struct 25 25 { 26 26 int spaceAlloc; // 栈占据的空间大小（最大大小） 27 27 int stackSize; // 栈的实际大小 28 28 tour_struct** list; 29 29 } stack_struct; 30 30 31 31 typedef struct 32 32 { 33 33 int spaceAlloc; // 队列占据的空间大小（最大大小） 34 34 int head; // 队头 35 35 int tail; // 队尾 36 36 bool full; // 是否队满 37 37 tour_struct** list; 38 38 } queue_struct; 39 39 40 40 const int INFINITY = 1000000; // 最大代价，用于初始化 tour 41 41 const int NO_CITY = -1; // 标记 tour->cityList 中尚未使用的部分 42 42 const int MAX_STRING = 1000; // 输出字符串最大长度 43 43 int n; // 城市总数 44 44 int nThread; // 线程数 45 45 city_t home_town = 0; // 起点和终点的编号 46 46 cost_t *digraph; // 读取的邻接表数据 47 47 tour_struct *best_tour; // 保存最优轨迹 48 48 queue_struct *queue; // 搜索队列 49 49 int queueSize; // 搜索队的大小 50 50 int initialTourCount; // 首次分配给各线程的搜索队列数量 51 51 tour_struct* allocTour(stack_struct* avail);// 必需的声明，不然该函数与栈操作函数混在一起 52 52 53 53 // 基本函数 54 54 void usage(char* prog_name)// 提示信息 55 55 { 56 56 fprintf(stderr, "\n<Usage> %s <nThread> <digraph file>\n", prog_name); 57 57 exit(0); 58 58 } 59 59 60 60 void readDigraph(FILE* digraph_file)// 从文件读取邻接表 61 61 { 62 62 int i, j; 63 63 fscanf_s(digraph_file, "%d", &n);// 第一行数字时城市数量 64 64 if (n <= 0) 65 65 { 66 66 fprintf(stderr, "\n<readDigraph> Number of city = %d\n", n); 67 67 exit(1); 68 68 } 69 69 digraph = (cost_t*)malloc(n * n * sizeof(cost_t)); 70 70 for (i = 0; i < n; i++) 71 71 { 72 72 for (j = 0; j < n; j++) 73 73 { 74 74 fscanf_s(digraph_file, "%d", &digraph[i * n + j]); 75 75 if (i == j && digraph[i * n + j] != 0) 76 76 { 77 77 fprintf(stderr, "\n<readDigraph> Diagonal element [%d, %d] = %d\n", i, j, digraph[i * n + j]); 78 78 exit(2); 79 79 } 80 80 else if (i != j && digraph[i * n + j] <= 0) 81 81 { 82 82 fprintf(stderr, "\n<readDigraph> Off-diagonal element [%d, %d] = %d\n", i, j, digraph[i * n + j]); 83 83 exit(3); 84 84 } 85 85 } 86 86 } 87 87 #ifdef DEBUG 88 88 printf("\n\t<readDigraph> Read map\n"); 89 89 #endif 90 90 } 91 91 92 92 void printDigraph(void)// 打印邻接表 93 93 { 94 94 int i, j; 95 95 printf("\n\t<printDigraph> Order = %d\nMatrix = \n", n); 96 96 for (i = 0; i < n; i++) 97 97 { 98 98 for (j = 0; j < n; j++) 99 99 printf("%2d ", digraph[i * n + j]); 100100 printf("\n"); 101101 } 102102 printf("\n"); 103103 } 104104 105105 void initeTour(tour_struct* tour, cost_t cost)// 初始化搜索轨迹，要求给出代价 106106 { 107107 int i; 108108 tour->cityList[0] = 0; 109109 for (i = 1; i <= n; tour->cityList[i++] = NO_CITY); 110110 tour->cost = cost; 111111 tour->count = 1; 112112 } 113113 114114 void copyTour(tour_struct* tour1, tour_struct* tour2)// 搜索轨迹 tour1 拷贝给 tour2 115115 { 116116 memcpy(tour2->cityList, tour1->cityList, (n + 1) * sizeof(city_t)); 117117 tour2->count = tour1->count; 118118 tour2->cost = tour1->cost; 119119 } 120120 121121 void printTour(int my_rank, tour_struct* tour, char* title)// 输出搜索轨迹，注意有两种输出 122122 { 123123 int i; 124124 char string[MAX_STRING]; 125125 if (my_rank >= 0)// 各线程汇报 126126 sprintf_s(string, "\n\t<printTour> Rank%2d, %s %p: ", my_rank, title, tour); 127127 else // 输出当前最佳回路 128128 sprintf_s(string, "\n\t<printTour> %s: ", title); 129129 for (i = 0; i < City_count(tour); i++) 130130 sprintf_s(string + strlen(string), 10, "%d ", Tour_city(tour, i)); 131131 printf("%s\n", string); 132132 } 133133 134134 // 栈相关 135135 stack_struct* initeStack(void)// 初始化栈 136136 { 137137 int i; 138138 stack_struct* stack = (stack_struct*)malloc(sizeof(stack_struct)); 139139 stack->list = (tour_struct**)malloc(n * n * sizeof(tour_struct*)); 140140 for (i = 0; i < n * n; stack->list[i++] = NULL); 141141 stack->stackSize = 0; 142142 stack->spaceAlloc = n * n; 143143 return stack; 144144 } 145145 146146 void pushStack(stack_struct* stack, tour_struct* tour)// 非拷贝压栈 147147 { 148148 if (stack->stackSize == stack->spaceAlloc) 149149 { 150150 fprintf(stderr, "\n<pushStack> Overflow\n"); 151151 free(tour->cityList); 152152 free(tour); 153153 } 154154 else 155155 { 156156 # ifdef DEBUG 157157 printf("\n\t<pushStack> StackSize = %d, tour at %p, tour->cityList at %p\n", stack->stackSize, tour, tour->cityList); 158158 printTour(-1, tour, "pushStack"); 159159 printf("\n"); 160160 # endif 161161 stack->list[stack->stackSize] = tour; 162162 (stack->stackSize)++; 163163 } 164164 } 165165 166166 void pushCopyStack(stack_struct* stack, tour_struct* tour, stack_struct* avail)// 拷贝压栈 167167 { 168168 tour_struct* tmp; 169169 if (stack->stackSize == stack->spaceAlloc) 170170 { 171171 fprintf(stderr, "\n<pushCopyStack> Overflow\n"); 172172 exit(-1); 173173 } 174174 tmp = allocTour(avail); 175175 copyTour(tour, tmp); 176176 stack->list[stack->stackSize] = tmp; 177177 (stack->stackSize)++; 178178 } 179179 180180 tour_struct* popStack(stack_struct* stack)// 出栈 181181 { 182182 tour_struct* tmp; 183183 if (stack->stackSize == 0) 184184 { 185185 fprintf(stderr, "\n<popStack> Empty\n"); 186186 exit(-1); 187187 } 188188 tmp = stack->list[stack->stackSize - 1]; 189189 stack->list[stack->stackSize - 1] = NULL; 190190 (stack->stackSize)--; 191191 return tmp; 192192 } 193193 194194 int isEmptyStack(stack_struct* stack)// 判断是否空栈 195195 { 196196 return stack->stackSize == 0; 197197 } 198198 199199 void freeStack(stack_struct* stack)// 清空栈 200200 { 201201 int i; 202202 for (i = 0; i < stack->stackSize; free(stack->list[i]->cityList), free(stack->list[i]), i++); 203203 free(stack->list); 204204 free(stack); 205205 } 206206 207207 void printStack(stack_struct* stack, int my_rank, char title[])// 打印栈 208208 { 209209 char string[MAX_STRING]; 210210 int i, j; 211211 printf("\n\t<printStack> Rank%2d, %s\n", my_rank, title); 212212 for (i = 0; i < stack->stackSize; i++) 213213 { 214214 sprintf_s(string, 10, "%d> ", i); 215215 for (j = 0; j < stack->list[i]->count; j++) 216216 sprintf_s(string + strlen(string), 10, "%d ", stack->list[i]->cityList[j]); 217217 printf("%s\n", string); 218218 } 219219 } 220220 221221 // 队列相关 222222 queue_struct* initeQueue(int size)// 初始化队列 223223 { 224224 queue_struct* new_queue = (queue_struct*)malloc(sizeof(queue_struct)); 225225 new_queue->list = (tour_struct**)malloc(size * sizeof(tour_struct*)); 226226 new_queue->spaceAlloc = size; 227227 new_queue->head = new_queue->tail = new_queue->full = false; 228228 return new_queue; 229229 } 230230 231231 int isEmptyQueue(queue_struct* queue)// 判断是否队空 232232 { 233233 return !queue->full && queue->head == queue->tail;// 空队要求 queue->full == false 且头尾指针相等 234234 } 235235 236236 tour_struct* deQueue(queue_struct* queue)// 出队 237237 { 238238 tour_struct* tmp; 239239 if (isEmptyQueue(queue)) 240240 { 241241 fprintf(stderr, "\n<deQueue> Empty queue\n"); 242242 exit(-1); 243243 } 244244 tmp = queue->list[queue->head]; 245245 queue->head = (queue->head + 1) % queue->spaceAlloc; 246246 return tmp; 247247 } 248248 249249 void enQueue(queue_struct* queue, tour_struct* tour)// 入队 250250 { 251251 tour_struct* tmp; 252252 if (queue->full == true) 253253 { 254254 fprintf(stderr, "\n<enQueue> Overflow\n"); 255255 exit(-1); 256256 } 257257 tmp = allocTour(NULL); 258258 copyTour(tour, tmp); 259259 queue->list[queue->tail] = tmp; 260260 queue->tail = (queue->tail + 1) % queue->spaceAlloc; 261261 if (queue->tail == queue->head) 262262 queue->full = true; 263263 } 264264 265265 void freeQueue(queue_struct* queue)// 清空队列 266266 { 267267 free(queue->list); 268268 free(queue); 269269 } 270270 271271 void printQueue(queue_struct* queue, int my_rank, char title[])// 打印队列 272272 { 273273 char string[MAX_STRING]; 274274 int i, j; 275275 printf("\n\t<printQueue> Rank%2d > %s\n", my_rank, title); 276276 for (i = queue->head; i != queue->tail; i = (i + 1) % queue->spaceAlloc) 277277 { 278278 sprintf_s(string, "%d> %p = ", i, queue->list[i]); 279279 for (j = 0; j < queue->list[i]->count; j++) 280280 sprintf_s(string + strlen(string), 10, "%d ", queue->list[i]->cityList[j]); 281281 printf("%s\n", string); 282282 } 283283 } 284284 285285 // Tour 相关 286286 tour_struct* allocTour(stack_struct* avail)// 生成一个搜索轨迹 287287 { 288288 tour_struct* tmp; 289289 if (avail == NULL || isEmptyStack(avail)) 290290 { 291291 tmp = (tour_struct*)malloc(sizeof(tour_struct)); 292292 tmp->cityList = (city_t*)malloc((n + 1) * sizeof(city_t)); 293293 return tmp; 294294 } 295295 return popStack(avail); 296296 } 297297 298298 int visited(tour_struct* tour, city_t city)// 判断一个轨迹中是否经过了城市 city 299299 { 300300 for (int i = 0; i < City_count(tour); i++) 301301 { 302302 if (Tour_city(tour, i) == city) 303303 return true; 304304 } 305305 return false; 306306 } 307307 308308 void addCity(tour_struct* tour, city_t new_city)// 将一个城市添加到 tour 309309 { 310310 city_t old_last_city = Last_city(tour); // 从原来的最后一个城市接出一条路径 311311 tour->cityList[tour->count] = new_city; 312312 (tour->count)++; 313313 tour->cost += Cost(old_last_city, new_city); 314314 } 315315 316316 void removeLastCity(tour_struct* tour)// 去掉 tour 中最后一个城市 317317 { 318318 city_t old_last_city = Last_city(tour), new_last_city; 319319 tour->cityList[tour->count - 1] = NO_CITY; 320320 (tour->count)--; 321321 new_last_city = Last_city(tour); 322322 tour->cost -= Cost(new_last_city, old_last_city); 323323 } 324324 325325 int isBestTour(tour_struct* tour)// 判断当前轨迹是否最佳 326326 { 327327 cost_t cost_so_far = Tour_cost(tour); 328328 city_t last_city = Last_city(tour); 329329 return cost_so_far + Cost(last_city, home_town) < Tour_cost(best_tour);// 当前轨迹的代价加上最后一个城市回家的代价是否更优 330330 } 331331 332332 int isFeasible(tour_struct* tour, city_t city)// 判断 tour 的最后一个城市是否可以前往城市 city 333333 { 334334 return !visited(tour, city) && Tour_cost(tour) + Cost(Last_city(tour), city) < Tour_cost(best_tour);// 要求没去过该城市且代价不超出当前最优回路代价 335335 } 336336 337337 void updateBestTour(tour_struct* tour) 338338 { 339339 if (isBestTour(tour))// 再次检查是否最优，防止两次检查之间其他线程更新了最优轨迹 340340 { 341341 copyTour(tour, best_tour); 342342 addCity(best_tour, home_town); 343343 } 344344 } 345345 346346 void freeTour(tour_struct* tour, stack_struct* avail) 347347 { 348348 if (avail == NULL) // 普通释放 349349 { 350350 free(tour->cityList); 351351 free(tour); 352352 } 353353 else // 将结点压回栈中 354354 pushStack(avail, tour); 355355 } 356356 357357 void distributeTour(int my_rank, int* my_first_tour_p, int* my_last_tour_p)// 分配初始搜索任务，使用块划分 358358 { 359359 const int quotient = initialTourCount / nThread, remainder = initialTourCount % nThread;// 平均每线程 quotient 个，余数部分由靠前的线程分担 360360 int my_count; 361361 if (my_rank < remainder)// 靠前的线程，分担余数 362362 { 363363 my_count = quotient + 1; 364364 *my_first_tour_p = my_rank*my_count; 365365 } 366366 else // 靠后的线程，直接分配，注意偏移 367367 { 368368 my_count = quotient; 369369 *my_first_tour_p = my_rank*my_count + remainder; 370370 } 371371 *my_last_tour_p = *my_first_tour_p + my_count - 1; 372372 } 373373 374374 // 核心计算函数附属 375375 inline long long factorial(int k)// 计算 k 的阶乘 376376 { 377377 long long tmp = 1; 378378 int i; 379379 for (i = 2; i <= k; tmp *= i++); 380380 return tmp; 381381 } 382382 383383 inline int supQueueSize(void)// 估计搜索需要的队列长度上限，并判断使用的线程数是否太多，没看懂 384384 { 385385 int fact, size; 386386 for (fact = size = n - 1; size < nThread; size *= ++fact); 387387 if (size > factorial(n - 1)) 388388 { 389389 fprintf(stderr, "\n<supQueueSize> Too many thread\n"); 390390 size = 0; 391391 } 392392 return size; 393393 } 394394 395395 void buildInitialQueue(void)// 生成初始搜索队列 396396 { 397397 int currentQueueSize; // 当前队列的长度 398398 city_t nbr; 399399 tour_struct* tour = allocTour(NULL); 400400 401401 initeTour(tour, 0); // 仅包含起点，代价为 0 402402 queue = initeQueue(2 * queueSize); // 搜索队 403403 404404 enQueue(queue, tour); // 将起点拷贝入队 405405 freeTour(tour, NULL); 406406 for (currentQueueSize = 1; queueSize < nThread;)// 向搜索队中添加元素，直到不小于线程数为止，以便分配各线程工作 407407 { 408408 tour = deQueue(queue); // 从队列中取出一个城市来 409409 currentQueueSize--; 410410 for (nbr = 1; nbr < n; nbr++) 411411 { 412412 if (!visited(tour, nbr)) // 取出的城市是新城市 413413 { 414414 addCity(tour, nbr); // 将该城市添加进 tour 415415 enQueue(queue, tour); // 将当前 tour 拷贝加入搜索队列中 416416 currentQueueSize++; 417417 removeLastCity(tour); // 加入搜索队之后立即从 tour 中清楚刚加入的城市，tour 仅用于搜索队操作 418418 } 419419 } 420420 freeTour(tour, NULL); // 本轮入队结束 421421 } 422422 initialTourCount = currentQueueSize;// 记录当前搜索队中 tour 的数量 423423 # ifdef DEBUG 424424 printQueue(queue, 0, "\n\t<buildInitialQueue> Queue Initialized\n"); 425425 # endif 426426 } 427427 428428 void treePartition(int my_rank, stack_struct* stack)// 分树 429429 { 430430 int my_first_tour, my_last_tour, i; 431431 # pragma omp single // 估计搜索队列长度上限 432432 queueSize = supQueueSize(); 433433 # ifdef DEBUG 434434 printf("\n\t<treePartition> Rank%2d> supQueueSize = %d\n", my_rank, queueSize); 435435 # endif 436436 if (queueSize == 0) // 线程数太多 437437 exit(0); 438438 439439 # pragma omp master // 主线程创建初始搜索队列 440440 buildInitialQueue(); 441441 # pragma omp barrier // 从线程要显式的等待主线程 442442 443443 distributeTour(my_rank, &my_first_tour, &my_last_tour);// 分配初始搜索任务 444444 # ifdef DEBUG 445445 printf("\n\t<treePartition> Rank%2d> initialTourCount = %d, first = %d, last = %d\n", my_rank, initialTourCount, my_first_tour, my_last_tour); 446446 # endif 447447 for (i = my_last_tour; i >= my_first_tour; i--)// 从分配到的任务列中中最后一个向前逐个压入栈中 448448 { 449449 # ifdef DEBUG 450450 printTour(my_rank, Queue_elt(queue, i), "pushStack"); 451451 # endif 452452 pushStack(stack, Queue_elt(queue, i)); 453453 } 454454 # ifdef DEBUG 455455 printStack(stack, my_rank, "Task set up"); 456456 # endif 457457 458458 } 459459 460460 // 核心计算函数 461461 void treeSearch(void) 462462 { 463463 int my_rank = omp_get_thread_num(); 464464 city_t nbr; 465465 stack_struct *stack = initeStack(), *avail = initeStack();// 搜索栈和未搜索的栈 466466 tour_struct *curr_tour; 467467 468468 for (treePartition(my_rank, stack); !isEmptyStack(stack);)// 初次分树，然后全部搜索完以前不断循环 469469 { 470470 curr_tour = popStack(stack); 471471 # ifdef DEBUG 472472 printTour(my_rank, curr_tour, "popStack"); 473473 # endif 474474 if (City_count(curr_tour) == n) // 当前回路已经搜索了 n 个城市 475475 { 476476 if (isBestTour(curr_tour)) // 判断是否是最佳回路 477477 { 478478 # ifdef DEBUG 479479 printTour(my_rank, curr_tour, "Best tour"); 480480 # endif 481481 # pragma omp critical 482482 updateBestTour(curr_tour);// 更新当前最佳回路 483483 } 484484 } 485485 else // 还没有搜索 n 个城市 486486 { 487487 for (nbr = n - 1; nbr >= 1; nbr--) 488488 { 489489 if (isFeasible(curr_tour, nbr)) 490490 { 491491 addCity(curr_tour, nbr); 492492 pushCopyStack(stack, curr_tour, avail); 493493 removeLastCity(curr_tour); 494494 } 495495 } 496496 } 497497 freeTour(curr_tour, avail); 498498 } 499499 freeStack(stack); // 释放线程栈资源 500500 freeStack(avail); 501501 # pragma omp barrier // 等待所有线程都完成 502502 # pragma omp master 503503 freeQueue(queue); // 主线程释放队列资源 504504 } 505505 506506 int main(int argc, char* argv[]) 507507 { 508508 FILE* digraph_file; // 输入文件名 509509 double start, finish; // 计时器 510510 # ifdef DEBUG 511511 argc = 3; 512512 argv[1] = "1"; 513513 argv[2] = "D:\\Code\\并行程序设计导论 - 代码\\ch6\\TSP\\mat_17e"; 514514 # endif 515515 516516 if (argc != 3) // 检查输入 517517 usage(argv[0]); 518518 if ((nThread = strtol(argv[1], NULL, 10)) <= 0) 519519 { 520520 fprintf(stderr, "\n<main> Error thread number\n"); 521521 usage(argv[0]); 522522 } 523523 fopen_s(&digraph_file, argv[2], "r"); 524524 if (digraph_file == NULL) 525525 { 526526 fprintf(stderr, "\n<main> Error opening input file\n"); 527527 usage(argv[0]); 528528 } 529529 530530 readDigraph(digraph_file);// 读取图 531531 fclose(digraph_file); 532532 # ifdef DEBUG 533533 printDigraph(); 534534 # endif 535535 536536 best_tour = allocTour(NULL); 537537 initeTour(best_tour, INFINITY); 538538 # ifdef DEBUG 539539 printTour(-1, best_tour, "Best tour"); 540540 printf("\n\t<main> City count = %d\nCost = %d\n\n", City_count(best_tour), Tour_cost(best_tour)); 541541 # endif 542542 543543 start = omp_get_wtime();// 开始计算 544544 # pragma omp parallel num_threads(nThread) default(none) 545545 treeSearch(); 546546 finish = omp_get_wtime(); 547547 548548 printTour(-1, best_tour, "Best tour");// 汇报结果 549549 printf("\n\t<main> Cost = %d\nElapsed time = %e s\n", best_tour->cost, finish - start); 550550 551551 free(best_tour->cityList);// 释放资源 552552 free(best_tour); 553553 free(digraph); 554554 return 0; 555555 }
● 输入的邻接表图
117 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 3 2 0 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 4 2 2 0 2 1 2 2 2 2 2 2 2 2 2 2 2 2 5 2 2 2 0 2 1 2 2 2 2 2 2 2 2 2 2 2 6 2 2 2 2 0 2 1 2 2 2 2 2 2 2 2 2 2 7 2 2 2 2 2 0 2 1 2 2 2 2 2 2 2 2 2 8 2 2 2 2 2 2 0 2 1 2 2 2 2 2 2 2 2 9 2 2 2 2 2 2 2 0 2 1 2 2 2 2 2 2 2 10 2 2 2 2 2 2 2 2 0 2 1 2 2 2 2 2 2 11 2 2 2 2 2 2 2 2 2 0 2 1 2 2 2 2 2 12 2 2 2 2 2 2 2 2 2 2 0 2 1 2 2 2 2 13 2 2 2 2 2 2 2 2 2 2 2 0 2 1 2 2 2 14 2 2 2 2 2 2 2 2 2 2 2 2 0 2 1 2 2 15 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 1 2 16 1 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 17 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 0 2 18 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0
● 输出结果，因为是指数级的枚举，没有等待程序运行结束