多線程動態規劃算法求解TSP(Traveling Salesman Problem) 並附C語言實現例程
TSP問題描述：
旅行商問題，即TSP問題（Travelling Salesman Problem）又譯為旅行推銷員問題、貨郎擔問題，是數學領域中著名問題之一。假設有一個旅行商人要拜訪n個城市，他必須選擇所要走的路徑，路徑的限制是每個城市只能拜訪一次，而且最後要回到原來出發的城市。路徑的選擇目標是要求得的路徑路程為所有路徑之中的最小值。這篇文章解決的tsp問題的輸入描述是：
TSP問題的動態規劃解法：
引用一下這篇文章，覺得作者把動態規劃算法講的非常明白：https://blog.csdn.net/derek_tian/article/details/45057443
動態規劃算法的並行：
終於到了文章的重點，如何高效率的將主問題拆開並分配給各個線程運算是值得研究的話題。這裏先引用動態規劃算法的示意圖：
（圖2）
圖中所示的是城市數為4的問題，從第一行到第二行，把問題拆分成了三個子問題，這三個子問題互不影響互不相關。將這三個子問題再一次拆分，每個問題又能拆出兩個問題。將這個規律推廣到n，城市數為n的問題通過第一次拆分可以變為n-1個子問題，再拆分一次共可以變為(n-1)*(n-2)個子問題。
有了能把問題拆成互不相幹的子問題作為前提，在並行時就並不需為數據頻繁的加鎖解鎖，能比較順利地完成並行的操作了。接下來是具體的實現手段：
首先在全局維護唯一一個二維表，所有子線程的讀取與修改都是在這一張表上進行的。例程中就是dp變量。
接下來定義一個函數，用來求解每個子問題，這個時候出現了另一個問題：怎麽樣表示一個子問題。在動態規劃算法中有一個狀態壓縮的概念，即為用一段二進制碼表示城市的集合，例如共五個城市時表示{1, 3, 4}這個集合即為二進制"11010"。相同的思路，對於d(1,{2,3})問題，我們給函數輸入起點，以及目標的集合，完整的表示問題。函數內選擇的是遞歸算法，因為非遞歸算法很難保證子問題互不相關（可能也只是我懶得想，有思路的同學可以在評論區中討論一下）。具體的動態規劃遞歸實現算法很多文章中都有，我就不過多贅述了。下面是這一部分的源代碼：
*註意：我在寫這個函數的時候規定的setMask是包含了起點城市的，也就是對於d(1,{2,3})這個問題傳入函數的參數的是(1,0b"1110")。
int TSP( int x, int setMask ) { int mask, i, minimum; if ( dp[ setMask ][ x ] != -1 ) return dp[ setMask ][ x ]; /* if already have value */ mask = 1 << x; dp[ setMask ][ x ] = 1e9; /* set infinite*/ for ( i = 0; i < n; ++i ) { if ( i != x && ( setMask & ( 1 << i ) ) && map[ i ][ x ] != 0 ) { /* if have path */ minimum = TSP( i, setMask - mask ) + map[ i ][ x ] ; dp[ setMask ][ x ] = min( dp[ setMask ][ x ], minimum); } } return dp[ setMask ][ x ]; }
現在我們有了表達子問題的方法，就該開始解決具體怎麽拆開問題的算法了。問題的要求中說明了最大可以創建的線程數，也就是說為了保證盡可能高的效率，程序中應該將問題拆到數量剛好大過線程數的情況。
當我們計算出了應該拆分的層數後，就該為產生的一大波子問題生成描述了。
有了之前定義子問題動態規劃的算法為基礎，其實可以發現，對於同一個setMask，其產生的同輩問題數是很容易得出來的。就拿剛剛的例子d(1,{2,3})來說，它的setMask是"1110"，起點是1號城市。然而，對於與這個問題互為同輩關系的d(2,{1,3})，d(3,{1,3})問題來說，他們的setMask也都是"1110"，只不過起點依次是2和3號城市。所以只要生成出了這個集合，其對應的問題都能很容易表示，這也就把工作的重心轉義到了如何生成這個集合，也就是例程中是setMask。
對於從同一個整道題求解的問題分離出來的子問題來說，其城市集合其實是一個組合數。我們還選取圖2作為例子，現在經過計算後，我們需要把主問題拆分兩層，即拆到d(2,{3})這一層。我們列出所有的集合，分別是"1100", "1010", "1010"，發現了嗎，其實這些集合也就是剔除了起點城市以後，將兩個1與一個0進行組合的所有情況。這個規律同樣可以推廣到n個城市，如果要將問題拆分k層，所有的問題集合即是n-k-1個1與k個0的所有組合數。
希望我講清楚了裏面的數學關系，如果覺得我講的有問題或者看不懂的同學可以在評論區留言……以下是上述數學過程的c語言代碼，為了邏輯更清晰我將計算組合數和求解所有組合數封成了factorial和setCombination函數：
while (t > nn) { /* let all thread have work to do */ nn *= (nn - 1); nc += 1; /* extends layers until t > nn */ } j = 0; candiNum = malloc(sizeof(int) * (n - 1)); for (i = 0; i < n; ++i) { /* candiNum set means all cities exclude the start piont*/ if (i != s) { candiNum[j] = i; /* exclude the start point */ j++; } } combNum = factorial(nc, n - 1); /* calculate combination number */ combSet = malloc(sizeof(int) * combNum); combptr = combSet; numOf1 = n - 1 - nc; /* there should have `numOf1` of 1 in every item */ setCombination(candiNum, n - 1, numOf1); posOf1 = malloc(sizeof(int) * numOf1); /* arr to store the position of 1 in each comb */
接下來就是使用pthread創建線程來求解每一個子問題了。其中主函數作為調度線程來給0~t-1個線程依次分配任務。因為基本上可以認為每一個平行的子問題計算的時間是完全一樣的，所以只需要先給0~t-1個線程分配0~t-1號任務，再等待到0號線程結束以後，為其分配t號任務，以此類推。當所有的子問題都被解決以後，dp這個二維表中已經有所有我們需要的子問題結論了。此時主問題的答案可以很快通過這些答案計算出來，也沒有並行的必要了。最後提供整個程序的代碼以及Makefile：
// file tsp.c
#include <stdio.h> #include <stdlib.h> #include <string.h> #include <pthread.h> #include "tsp.h" #include <time.h> #define min(a,b) (((a)<(b))?(a):(b)) /* macro func of min */ time_t tm; int ** map, ** dp, * combSet, * combptr, * posOf1, n; /* init vars */ trdArgs args[20]; pthread_t tid[20]; /* thread ids */ int TSP( int x, int setMask ) { int mask, i, minimum; if ( dp[ setMask ][ x ] != -1 ) return dp[ setMask ][ x ]; /* if already have value */ mask = 1 << x; dp[ setMask ][ x ] = 1e9; /* set infinite */ for ( i = 0; i < n; ++i ) { if ( i != x && ( setMask & ( 1 << i ) ) && map[ i ][ x ] != 0 ) { /* if have path */ minimum = TSP( i, setMask - mask ) + map[ i ][ x ] ; dp[ setMask ][ x ] = min( dp[ setMask ][ x ], minimum); } } return dp[ setMask ][ x ]; } void * trd(void * SArg) { TSP((*(trdArgs *)SArg).startPoint, (*(trdArgs *)SArg).setMask); return NULL; } void setCombination(int arr[], int n, int r) { /* set all combination into combptr */ int * data = malloc(sizeof(int) * r); /* Temporary array to store current combination */ combinationUtil(arr, data, 0, n - 1, 0, r); free(data); /* avoid leak */ } void combinationUtil(int arr[], int data[], int start, int end, int index, int r) { int i, j, tmpi = 0; /* init vars */ if (index == r) { for (j = 0; j < r; j++) /* add nums to binay */ tmpi |= (1 << data[j]); *combptr = tmpi; combptr++; /* next cell in array */ return; } for (i = start; i <= end && end - i + 1 >= r - index; i++) { /* recersive */ data[index] = arr[i]; combinationUtil(arr, data, i + 1, end, index + 1, r); /* call selfe */ } } long factorial(int m, int n) { /* caculate combination number */ int i, j; long ans = 1; if (m < n - m) m = n - m; /* C(m,n)=C(n-m,n) */ for (i = m + 1; i <= n; i++) ans *= i; for (j = 1; j <= n - m; j++) ans /= j; return ans; /* answer */ } int main() { long combNum; int t, s, i, j, k, nn, nc, ans, numOf1, * candiNum, * combSet, * tmptr; /* init vars */ nc = 0; scanf("%d %d %d", &t, &n, &s); /* get &t, &n, &s */ if (t < 0 || n <= 0 || n <= s || s < 0) {printf("-1\n"); return 0;} /* error input */ if (n == 1) {printf("0\n"); return 0;} nn = n - 1; /* first layer */ map = malloc(sizeof(int *)*n); /* alloc for 2dim array: map */ tmptr = malloc(sizeof(int) * n * n); for (i = 0; i < n; ++i) { map[i] = tmptr; /* set 2dim to 1dim */ tmptr += n; } dp = malloc(sizeof(int *) * (1 << n)); /* alloc for 2dim array: map */ tmptr = malloc(sizeof(int) * (1 << n) * n); for (i = 0; i < (1 << n); ++i) { dp[i] = tmptr; /* set 2dim to 1dim */ tmptr += n; } for (i = 0; i < n; i++) /* get the map */ for (j = 0; j < n; j++) scanf("%d", &map[i][j]); /* store */ memset(dp[0], -1, sizeof(int) * (1 << n) * n); /* fill all dp with -1 */ for ( i = 0; i < n; ++i ) dp[ 1 << i ][ i ] = map[ 0 ][ i ]; while (t > nn) { /* let all thread have work to do */ nn *= (nn - 1); nc += 1; /* extends layers until t > nn */ } j = 0; candiNum = malloc(sizeof(int) * (n - 1)); for (i = 0; i < n; ++i) { /* candiNum set means all cities exclude the start piont*/ if (i != s) { candiNum[j] = i; /* exclude the start point */ j++; } } combNum = factorial(nc, n - 1); /* calculate combination number */ combSet = malloc(sizeof(int) * combNum); combptr = combSet; numOf1 = n - 1 - nc; /* there should have `numOf1` of 1 in every item */ setCombination(candiNum, n - 1, numOf1); posOf1 = malloc(sizeof(int) * numOf1); /* arr to store the position of 1 in each comb */ k = 0; for (i = 0; i < combNum; ++i) { tmptr = posOf1; for (j = 0; j < n; ++j) { /* go through all bits */ if ((combSet[i] & (1 << j)) != 0) { *tmptr = j; tmptr++; /* point to next cell */ } } for (j = 0; j < numOf1; ++j) { /* arrange jobs */ args[k].id = k; /* set thread id arg */ args[k].startPoint = posOf1[j]; args[k].setMask = combSet[i]; /* set thread set mask args */ pthread_join(tid[k], NULL); pthread_create(&tid[k], NULL, trd, &args[k]); /* new thread */ k = (k + 1) % t; } } for (i = 0; i < t; ++i)/* join all thread */ pthread_join(tid[i], NULL); ans = TSP(0, ( 1 << n ) - 1); if (ans == 1e9) printf("-1\n"); else printf("%d\n", ans); /* final answer */ free(dp[0]); /* free */ free(dp); free(map[0]); /* free */ free(map); free(posOf1); /* free */ free(combSet); free(candiNum); /* free */ return 0; }
// file tsp.h typedef struct threadArgs { int id; int startPoint; int setMask; } trdArgs; int TSP( int x, int setMask ); void * trd(void * SArg) ; void combinationUtil(int arr[], int data[], int start, int end, int index, int r); void setCombination(int arr[], int n, int r) ; long factorial(int m, int n) ;
Makefile：
CC=gcc CFLAGS=-Wpedantic -Wall -Wextra -Werror -std=c89 -pthread -D_GNU_SOURCE tsp.o: tsp.c tsp.h ${CC} ${CFLAGS} tsp_blog.c -o tsp.o
歡迎沒有完全看懂的讀者隨時提問，文中的例程也不一定是最好的算法，還請有想法的讀者能不吝賜教。
多線程動態規劃算法求解TSP(Traveling Salesman Problem) 並附C語言實現例程