Parellel TSP

学校给了一个4×3G的cluster，写一个并行的程序。这里采用的是MPI， Message Passing Interface。算法在前面的post里讨论过了，我的算法很烂拉，没有好好利用，就算在那么强悍的机器上也没有办法跑到什么好结果。班上有4组牛人跑到最好的4个结果，166106，167692，171233，171313。哪位大哥有兴趣可以挑战一下，嘿嘿。

这里是代码，non-blocking的版本。写的很烂

tsp.h

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <mpi.h>
#include <time.h>

#ifndef _TSP_H_
#define _TSP_H_

#define CITYSIZE 500
/* 15000 > distance(0,0,10000,10000) */
#define MAXDIST  15000
/* 7500000 > 15000 * (500 - 1) */
#define MAXTOTAL 7500000

/* cities informatin */
typedef struct S_c {
 int x;
 int y;
 int visited;
} CITY;

/* a linked-list record local path generated by perm() */
struct S_path {
 int* path;
 struct S_path* next;
 struct S_path* prev;
} *head, *tail;

//***********************************************************************
void showCity(CITY*,int);              // show city coordinate,for test *
void showDist(int[][CITYSIZE],int);          // show distances,for test *
void showPath(int*,int,int,int[][CITYSIZE]);         // show found path *
double distance(int,int,int,int);                   // distance formula *
void perm(int*,int,int,int);                  // generate a permutation *
void swap(int*,int*);                                     // for perm() *
int localSearch(int,int*,int*,int[][CITYSIZE]);                      // *
int findrepeat(int*,int);                      // check if result valid *
int weightOf(int*,int[][CITYSIZE]);        // calculate a path's weight *
void ga(int);                               // simple genetic algorithm *
/************************************************************************/

void showDist(int array[][CITYSIZE], int size){
 int i, j;
 for(i = 0; i < size; ++i){
  for(j = 0; j < size; ++j){
   printf("%4d,", array[i][j]);
  }
  printf("/n");
 }
}

void showCity(CITY* array, int size){
 int i; 
 for(i = 0; i < size; ++i){
  printf("%d:%d:%d:%d/n", i,
        array[i].x,
        array[i].y,
        array[i].visited);
 }// for end
}

void showPath(int* array, int size, int total, int weight[][CITYSIZE]){
 int i;
 if(findrepeat(array, size)){printf("Wrong answer!/n");exit(0);}
 if(total!=weightOf(array, weight)){
  printf("wrong answer!/n");exit(0);}
 printf("%d:", total);
 for(i = 0; i < size - 1; ++i) printf("%d,", array[i]);
 printf("%d/n/n", array[i]);
}

double distance(int x1, int y1, int x2, int y2){
 return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

int localSearch(int step, int* tour, int* total, int weight[][CITYSIZE]){
 int i = 0;
 
 for(; (i + step) < CITYSIZE; i = i + step){
  int j = 0;
  int sublength = 0, newweight = 0, currBest;
  for(; j < step; ++j) sublength += weight[tour[i+j]][tour[i+j+1]];
  currBest = sublength;
  /* for sequence 01234, we need a permutation of p(123) only.
   * so, the lower and upper bound would be i + 1 and i + step - 1.
   */
  perm(tour, i + 1, i + step - 1, i + 1);
  // check all new path in permutation.
  struct S_path* onePath = tail->prev;
  while(onePath != head){
   onePath->prev->next = tail;
   tail->prev = onePath->prev;
   newweight = weight[tour[i]][onePath->path[0]] +
               weight[onePath->path[step-2]][tour[i+step]];
   for(j = 0; j < step - 2; ++j)
    newweight += weight[onePath->path[j]][onePath->path[j+1]];
   if(newweight < currBest){
    currBest = newweight;
    memcpy(tour+i+1,onePath->path,sizeof(int)*(step-1));
   }
   free(onePath->path);
   free(onePath);
   onePath = tail->prev;
  }// while end
  // update total length
  *total = *total - sublength + currBest;
 }// for end
 //showPath(tour,CITYSIZE, *total, weight);
 return *total;
}

/* variable 'original' record the start point(lower bound) in array,
 * used for output, when invoke perm(), original and k is same.
 */
void perm(int* list, int k, int m, int original){
 int i, j;
 if(k > m){// exist a permutation, add to list
  struct S_path* myPath =(struct S_path*)malloc(sizeof(struct S_path));
  struct S_path* tmp = head->next;
  myPath->path = (int*)malloc(sizeof(int)*(m - original + 1));
  head->next = myPath;
  myPath->prev = head;
  myPath->next = tmp;
  tmp->prev = myPath;
  for(i = original; i <= m; ++i) myPath->path[i-original] = list[i];
 } else {
  for(j = k; j <= m; ++j){
   swap(&list[k], &list[j]);
   perm(list, k + 1, m, original);
   swap(&list[k], &list[j]);
  }
 }
}

void swap(int* a, int* b){
 int temp = *a;
 *a = *b;
 *b = temp;
}

int findrepeat(int* array, int size){
 int i, j;
 for(i = 0; i < size; ++i){
  int tmp = 0;
  for(j = 0; j < size; ++j) if(i == array[j]) tmp++;
  if(tmp > 1){printf("error: repeat city!/n");return 1;}
 }
 return 0;
}

int weightOf(int* array, int weight[][CITYSIZE]){
 int w = 0, i = 0;
 for(; i < CITYSIZE - 1; ++i)
  w += weight[array[i]][array[i+1]];
 return w;
}

/*
// cannot get a path shorter than 190k in 1min, DEPRECATED
// need weight[][],all500paths[][],all500weights[] support
void ga(int separator){
 time_t seed;
 time(&seed);
 srand(seed);

 int epoch = 0;
 for(; epoch < 1000000; ++epoch){
  // printf("epoch %d:/n", epoch);
  rand(); // dump first number
  
  // randomly pick two path index
  int index1 = (int)((double)rand()*500/RAND_MAX);
  int index2 = (int)((double)rand()*500/RAND_MAX);

  // get the two paths copies
  int path1[CITYSIZE], path2[CITYSIZE];
  memcpy(path1, all500paths[index1], sizeof(path1));
  memcpy(path2, all500paths[index2], sizeof(path2));

  // exchange two fragments to get next epoch
  int* frag1 = (int*)malloc(sizeof(int)*separator);
  memcpy(frag1, path1, sizeof(int)*separator);
  memcpy(path1, path2, sizeof(int)*separator);
  memcpy(path2, frag1, sizeof(int)*separator);
  free(frag1);

  // re-generate path1 & path2
  int i, check1[CITYSIZE] = {0}, check2[CITYSIZE] = {0};
  // separator~CITYSIZE-1 is original path, isn't changed
  for(i = separator; i < CITYSIZE; ++i){
   check1[path1[i]]++;
   check2[path2[i]]++;
  }

  // check new part in path
  // 0~seprator-1 is new part copied from the other one
  for(i = 0; i < separator; ++i){
   if(check1[path1[i]]){
    int j = 0;
    while(check1[j])j++; // step to first non-zero position
    check1[j]++;
    path1[i] = j;
   }else check1[path1[i]]++;
  
   if(check2[path2[i]]){
    int j = 0;
    while(check2[j])j++;
    check2[j]++;
    path2[i] = j;
   }else check2[path2[i]]++;
  }

  int mytotal1 = 0, mytotal2 = 0;
  // calculate new path length
  for(i = 0; i < CITYSIZE - 1; ++i){
   mytotal1 += weight[path1[i]][path1[i+1]];
   mytotal2 += weight[path2[i]][path2[i+1]];
  }
  // update if better route found
  if(mytotal1 < all500weights[path1[0]]){
   printf("better found in path%d: %d/n", path1[0], mytotal1);
   all500weights[path1[0]] = mytotal1;
   memcpy(all500paths[path1[0]], path1, sizeof(int)*CITYSIZE);
  }
  if(mytotal2 < all500weights[path2[0]]){
   printf("better found in path%d: %d/n", path2[0], mytotal2);
   all500weights[path2[0]] = mytotal2;
   memcpy(all500paths[path2[0]], path2, sizeof(int)*CITYSIZE);
  }
 }// for end, one epoch end
}
*/
#endif

tsp.c

/* 500 cities Hamilton Path */
/* Repetitive Nearest Neighbour Algorithm (RNNA) */
/* with Local Search */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <mpi.h>
#include <time.h>
#include "tsp.h"

/* local search factor
 * greater number means deeper recursive
 */
#define STEP 11

#define WORK 1
#define DIE 2

int tour[CITYSIZE],  // best path
 total = MAXTOTAL; // best path length

void master();
void slave(int);

int main(int argc, char *argv[]){
// time_t start, done;
// time(&start);
    int myid;

 MPI_Init(&argc, &argv);
 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
   
 if(myid == 0) master();
 else slave(myid);

    MPI_Finalize();
// time(&done);
// printf("process %d's execution time = %d seconds/n", myid, done - start);
    return 0;
}

void master(){
 // mpi used varibles
 int numprocs;
 int allFree = 1;
 char tmp;
 int* freeprocs = (int*)calloc(numprocs - 1, sizeof(int));
 int* indicies  = (int*)malloc((numprocs - 1) * sizeof(int));
 MPI_Request* reqs = (MPI_Request*)calloc(numprocs-1, sizeof(MPI_Request));
 MPI_Status * stat = (MPI_Status*)malloc(sizeof(MPI_Status)*(numprocs-1));
 MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
 memset(freeprocs, 0, numprocs - 1);
 
 // assignment used varibles
 FILE *fin;
 CITY city[CITYSIZE];
 int i, j,
  tmptour[CITYSIZE],
  weight[CITYSIZE][CITYSIZE]={0};   // distances between cities

 fin = fopen("cities.dat", "r");
 
 /* read cities locations */
 char dump;  // useless information, used to omit ':'
 for(i = 0; i < CITYSIZE; ++i){
  fscanf(fin, "%d %c %d", &dump, &dump, &city[i].x);
  fscanf(fin, "%c %d", &dump, &city[i].y);
 }
 
 /* calculate cities distance */
 for(i = 0; i < CITYSIZE; ++i){
     for(j = 0; j < CITYSIZE; ++j){
      if(i < j){
          weight[i][j] = (int)distance(city[i].x, city[i].y,
                                    city[j].x, city[j].y);
          weight[j][i] = weight[i][j];
      }//if end
     }//for end
 }// for end

 /* RNNA start, use NNA for all 500 cities.
  * best path stored in tour[], best weight stored in total
  */
 int tmptotal = MAXTOTAL, startCty = 0;
 for(; startCty < CITYSIZE; ++startCty){
  /* if found better path */
  if(tmptotal < total){
   total = tmptotal;
   memcpy(tour, tmptour, sizeof(tour));
   showPath(tour, CITYSIZE, total, weight);
  }
  /* initialization */
  tmptotal = 0;
  for(i = 0; i < CITYSIZE; ++i) city[i].visited = 0;
  
  int nextCty, currCty = startCty;
  tmptour[0] = startCty;
  city[startCty].visited = 1;
  
  for(i = 1; i < CITYSIZE; ++i){
   int cost = MAXDIST;
   for(j = 0; j < CITYSIZE; ++j){
    if((cost > weight[currCty][j]) && // currCty~j is better and
       (!city[j].visited)  &&         // j hasn't been visited and
       (j != currCty)){               // j is not currCty
        cost = weight[currCty][j];
        nextCty = j;
    }// if end
   }// for end

   // step to next city
   currCty = nextCty;
   // add this city into path
   tmptour[i] = currCty;
   city[currCty].visited = 1;
   // accumulate path length
   tmptotal += cost;
  }// for end
 }// for end 
 /* RNNA finished here, local search start */
 
 // send distance matrix to slaves first
 for(i = 1; i < numprocs; ++i){
  MPI_Isend(weight, CITYSIZE * CITYSIZE, MPI_INT,
      i, WORK, MPI_COMM_WORLD, reqs + i);
 }
 
 // get packet length
 int s1, s2, s3, st;
 MPI_Pack_size(CITYSIZE, MPI_INT, MPI_COMM_WORLD, &s1);
 MPI_Pack_size(1, MPI_INT, MPI_COMM_WORLD, &s2);
 MPI_Pack_size(1, MPI_INT, MPI_COMM_WORLD, &s3);

 st = s1 + s2 + s3;
 
 int step = 3, position = 0;
 // used to store received paths
 int lsResults[numprocs - 1][CITYSIZE];
 while(step < STEP || allFree){
  for(i = 0; i < numprocs - 1 && step < STEP; ++i){
   if(*(freeprocs + i) == 0){// found idle slave...
    *(freeprocs + i) = 1;
    char* buffer = (char*)malloc(sizeof(char)*(st));
    MPI_Pack(tour, CITYSIZE, MPI_INT, buffer,
       st, &position, MPI_COMM_WORLD);
    MPI_Pack(&total, 1, MPI_INT, buffer,
       st, &position, MPI_COMM_WORLD);
    MPI_Pack(&step, 1, MPI_INT, buffer,
       st, &position, MPI_COMM_WORLD);
  
    MPI_Send(buffer, position, MPI_PACKED,
        i + 1, WORK, MPI_COMM_WORLD);
    free(buffer);
    position =0;
    step++;
    MPI_Irecv(*(lsResults + i), CITYSIZE, MPI_INT,
        i + 1, 1, MPI_COMM_WORLD, reqs + i);
   }// if end
  }// for end
  MPI_Testsome(numprocs - 1, reqs, &j, indicies, stat);

  if(j > 0){
   for(i = 0; i < j; ++i){
    *(freeprocs + indicies[i]) = 0;
    tmptotal = weightOf(lsResults[indicies[i]], weight);
    if(tmptotal < total){
     total = tmptotal;
     memcpy(tour, lsResults[indicies[i]], sizeof(tour));
     showPath(tour, CITYSIZE, total, weight);
    }// if end
   }// for end
  }// if end
  
  allFree = 1;
  for(i = 0; i < numprocs-1 && allFree; ++i)
   if(*(freeprocs + i) == 0) allFree = 0;
 }// while end
 
 // kill all slaves
 for(i = 1; i < numprocs; ++i){
  // no meaning data, just for matching data format
  char* buffer = (char*)malloc(sizeof(char)*(st));
  MPI_Pack(tour, CITYSIZE, MPI_INT, buffer,
     st, &position, MPI_COMM_WORLD);
  MPI_Pack(&total, 1, MPI_INT, buffer,
     st, &position, MPI_COMM_WORLD);
  MPI_Pack(&step, 1, MPI_INT, buffer,
     st, &position, MPI_COMM_WORLD);
  MPI_Isend(buffer, position, MPI_PACKED,
      i, DIE, MPI_COMM_WORLD, reqs + i);
  free(buffer);
  position = 0;
 }
}

void slave(int myid){
 MPI_Status status;
 MPI_Request request;
 int W[CITYSIZE * CITYSIZE],
  weight[CITYSIZE][CITYSIZE],
  tour[CITYSIZE],
  i, j, total, step, position = 0;

 // get packet length
 int s1, s2, s3, st;
 MPI_Pack_size(CITYSIZE, MPI_INT, MPI_COMM_WORLD, &s1);
 MPI_Pack_size(1, MPI_INT, MPI_COMM_WORLD, &s2);
 MPI_Pack_size(1, MPI_INT, MPI_COMM_WORLD, &s3);

 st = s1 + s2 + s3;
 
 // initialize list for localSearch()
 head = (struct S_path*)malloc(sizeof(struct S_path));
 tail = (struct S_path*)malloc(sizeof(struct S_path));
 head->next = tail;
 head->prev = NULL;
 tail->prev = head;
 tail->next = NULL;

 // receive distance matrix
 MPI_Recv(W, CITYSIZE * CITYSIZE, MPI_INT,
    0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
 // cast to 2d array
 for(i = 0; i < CITYSIZE; ++i)
  for(j = 0; j < CITYSIZE; ++j)
   weight[i][j] = W[i * CITYSIZE + j];
 
 while(1){
  char* buffer = (char*)malloc(sizeof(char)*st);
  MPI_Recv(buffer, st, MPI_PACKED,
     0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
  MPI_Unpack(buffer, st, &position,
       tour, CITYSIZE, MPI_INT, MPI_COMM_WORLD); 
  MPI_Unpack(buffer, st, &position,
       &total, 1, MPI_INT, MPI_COMM_WORLD);
  MPI_Unpack(buffer, st, &position,
       &step, 1, MPI_INT, MPI_COMM_WORLD);
  free(buffer);
  position = 0;
  if (status.MPI_TAG == DIE) return;
  
  localSearch(step, tour, &total, weight);
  MPI_Isend(tour, CITYSIZE, MPI_INT, 0,
      WORK, MPI_COMM_WORLD, &request);
  MPI_Wait(&request, &status);
 }
}