/* File:     pth_tsp_static.c
 * Purpose:  Implement a parallel tree search solution
 *           to an instance of the travelling salesman problem.
 *           This version uses static partitioning:  if there
 *           are n cities and p threads, thread 0 will search
 *           paths starting with 0->1, 0->2, . . . , 0->(n-1)/p,
 *           thread  1 will search paths starting with
 *           0->n/p+1, . . . , 0->2(n-1)/p, etc.
 *
 * Input:    From a user-specified file, the number of cities
 *           followed by the costs of travelling between the
 *           cities organized as a matrix:  the cost of
 *           travelling from city i to city j is the ij entry.
 * Output:   The best tour found by the program and the cost
 *           of the tour.
 *
 * Compile:  gcc -g -Wall -o pth_tsp pth_tsp_static.c -lpthread
 * Usage:    pth_tsp <thread count> <matrix_file>
 *
 * Notes:
 * 1.  The number of cities should be greater than the number
 *     of threads.  If it isn't, the program reduces the
 *     number of threads to (number of cities)-1.
 * 2.  Weights and cities are non-negative ints.
 * 3.  Program assumes the cost of travelling from a city to
 *     itself is zero, and the cost of travelling from one
 *     city to another city is positive.
 * 4.  The stacks are implemented with dynamic arrays.
 */
#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>

#define INFINITY 10000
#define NO_CITY -1
#define FALSE 0
#define TRUE 1
#define INITIAL_STACK_SIZE 1000

typedef int city_t;
typedef int weight_t;

typedef struct {
   city_t* cities;
   int count;
   weight_t cost;
} tour_t;

typedef struct {
   tour_t* tour_p;
   city_t  city;
} stack_elt_t;

typedef struct {
   stack_elt_t**   list;
   int             top;
   int             allocated;
/* pthread_mutex_t mutex;  Not needed in this version */
} my_stack_t;


/*------------------------------------------------------------------*/
/* Global variables */
int             thread_count;
int             n;
weight_t*       mat;
tour_t          best_tour;
pthread_mutex_t best_tour_mutex;
/* my_stack_t**       stack_ps; */

/*------------------------------------------------------------------*/
void Usage(char* prog_name);
void Read_mat(FILE* mat_file);
void Print_mat(void);

void Initialize_tour(tour_t* tour_p);
void* Search(void* rank);
void Print_tour(tour_t* tour_p, char* title);
void Check_best_tour(city_t city, tour_t* tour_p,
       weight_t* loc_best_cost);
int  Feasible(city_t city, city_t nbr, tour_t* tour_p,
      weight_t loc_best_cost);
int  Visited(city_t nbr, tour_t* tour_p);
void Initialize_stack(my_stack_t* stack_p);
void Push(tour_t* tour_p, city_t city, my_stack_t* stack_p);
tour_t* Dup_tour(tour_t* tour_p);
void Pop(tour_t** tour_pp, city_t* city_p, my_stack_t* stack_p);
int  Empty(my_stack_t* stack_p);
void Free_stack(my_stack_t* stack_p);
void Create_initial_records(long my_rank, my_stack_t* stack_p);
void Find_initial_cities(long my_rank, int* my_first_city_p, 
   int* my_last_city_p);


/*------------------------------------------------------------------*/
int main(int argc, char* argv[]) {
   long        thread;
   FILE*       mat_file;
   pthread_t*  thread_handles;

   if (argc != 3) Usage(argv[0]);
   thread_count = strtol(argv[1], NULL, 10);
   mat_file = fopen(argv[2], "r");
   if (mat_file == NULL) {
      fprintf(stderr, "Can't open %s\n", argv[1]);
      Usage(argv[0]);
   }
   Read_mat(mat_file);
   fclose(mat_file);
#  ifdef DEBUG
   Print_mat();
#  endif   
   if (thread_count >= n) {
      printf("Thread count should be < number of cities\n");
      printf("Changing thread count to %d\n", n-1);
      thread_count = n-1;
   }

   thread_handles = malloc(thread_count*sizeof(pthread_t));
/*
 * stack_ps = (my_stack_t**) malloc(thread_count*sizeof(my_stack_t*));
 * for (thread = 0; thread < thread_count; thread++)
 *    stack_ps[thread] = NULL;
 */
   pthread_mutex_init(&best_tour_mutex, NULL);
   Initialize_tour(&best_tour);
   best_tour.cost = INFINITY;

   for (thread = 0; thread < thread_count; thread++)
      pthread_create(&thread_handles[thread], NULL, Search,
         (void*) thread);
   for (thread = 0; thread < thread_count; thread++)
      pthread_join(thread_handles[thread], NULL);

   Print_tour(&best_tour, "Best tour");
   printf("Cost = %d\n", best_tour.cost);

   pthread_mutex_destroy(&best_tour_mutex);
   free(thread_handles);
   free(best_tour.cities);
   free(mat);
   return 0;
}  /* main */

/*------------------------------------------------------------------
 * Function:  Usage
 * Purpose:   Inform user how to start program and exit
 * In arg:    prog_name
 */
void Usage(char* prog_name) {
   fprintf(stderr, "usage: %s <thread_count> <matrix file>\n", prog_name);
   exit(0);
}  /* Usage */

/*------------------------------------------------------------------
 * Function:         Read_mat
 * Purpose:          Read in the number of cities and the matrix of costs
 * In arg:           mat_file
 * Global vars out:  mat, n
 */
void Read_mat(FILE* mat_file) {
   int i, j;

   fscanf(mat_file, "%d", &n);
   mat = malloc(n*n*sizeof(weight_t));

   for (i = 0; i < n; i++)
      for (j = 0; j < n; j++)
         fscanf(mat_file, "%d", &mat[n*i+j]);
}  /* Read_mat */


/*------------------------------------------------------------------
 * Function:        Print_mat
 * Purpose:         Print the number of cities and the matrix of costs
 * Global vars in:  mat, n
 */
void Print_mat(void) {
   int i, j;

   printf("Order = %d\n", n);
   printf("Matrix = \n");
   for (i = 0; i < n; i++) {
      for (j = 0; j < n; j++)
         printf("%2d ", mat[i*n+j]);
      printf("\n");
   }
   printf("\n");
}  /* Print_mat */


/*------------------------------------------------------------------*/
void Initialize_tour(tour_t* tour_p) {
   int i;

   tour_p->cities = malloc((n+1)*sizeof(city_t));
   for (i = 0; i <= n; i++) {
      tour_p->cities[i] = NO_CITY;
   }
   tour_p->cost = 0;
   tour_p->count = 0;
}  /* Initialize_tour */


/*------------------------------------------------------------------
 * Function:            Search
 * Purpose:             Search for an optimal tour
 * Global vars in:      mat, n
 * Global vars in/out:  best_tour
 */
void* Search(void* rank) {
   long     my_rank = (long) rank;
   city_t   nbr;
   city_t   city;
   tour_t*  tour_p;
   my_stack_t  stack;
   weight_t loc_best_cost = INFINITY;

   Initialize_stack(&stack);
   Create_initial_records(my_rank, &stack);
/* stack_ps[my_rank] = &stack; */

   while (!Empty(&stack)) {
      Pop(&tour_p, &city, &stack);
      tour_p->cities[tour_p->count] = city;
      tour_p->count++;
      if (tour_p->count == n) {
         Check_best_tour(city, tour_p, &loc_best_cost);
      } else {
         for (nbr = 1; nbr < n; nbr++) 
            if (Feasible(city, nbr, tour_p, loc_best_cost)) {
               tour_p->cost += mat[n*city+nbr];
               Push(tour_p, nbr, &stack);
               tour_p->cost -= mat[n*city+nbr];
            }
      }
      free(tour_p->cities);
      free(tour_p);
   }  /* while */

   Free_stack(&stack);

   return NULL;
}  /* Search */


/*------------------------------------------------------------------
 * Function:        Feasible
 * Purpose:         Check whether nbr could possibly lead to a better
 *                  solution if it is added to the current tour.  The
 *                  functions checks whether nbr has already been visited
 *                  in the current tour, and, if not, whether adding the
 *                  edge from the current city to nbr will result in
 *                  a cost less than the current best cost.
 * In args:         All
 * Global vars in:  mat, n, best_tour
 * Return:          TRUE if the nbr can be added to the current tour.
 *                  FALSE otherwise
 */
int Feasible(city_t city, city_t nbr, tour_t* tour_p, 
      weight_t loc_best_cost) {
   if (!Visited(nbr, tour_p) &&
        tour_p -> cost + mat[n*city + nbr] < loc_best_cost) 
      return TRUE;
   else
      return FALSE;
}  /* Feasible */


/*------------------------------------------------------------------
 * Function:   Visited
 * Purpose:    Use linear search to determine whether nbr has already
 *             bee visited on the current tour.
 * In args:    All
 * Return val: TRUE if nbr has already been visited.  
 *             FALSE otherwise
 */
int Visited(city_t nbr, tour_t* tour_p) {
   int i;

   for (i = 0; i < tour_p->count; i++)
      if ( tour_p->cities[i] == nbr ) return TRUE;
   return FALSE;
}  /* Visited */


/*------------------------------------------------------------------
 * Function:  Print_tour
 * Purpose:   Print a tour
 * In args:   All
 */
void Print_tour(tour_t* tour_p, char* title) {
   int i;

   printf("%s:\n", title);
   for (i = 0; i < tour_p->count; i++)
      printf("%d ", tour_p->cities[i]);
   printf("\n\n");
}  /* Print_tour */


/*------------------------------------------------------------------
 * Function:            Check_best_tour
 * Purpose:             Determine whether the current n-city tour will be
 *                      better than the current best tour.  If so, update
 *                      best_tour
 * In args:             city, tour_p
 * Global vars in:      mat, n
 * Global vars in/out:  best_tour
 */
void Check_best_tour(city_t city, tour_t* tour_p, 
       weight_t* loc_best_cost_p) {
   int i;

   pthread_mutex_lock(&best_tour_mutex);
   if (tour_p->cost + mat[city*n + 0] < best_tour.cost) {
      for (i = 0; i < tour_p->count; i++)
         best_tour.cities[i] = tour_p->cities[i];
      best_tour.cities[n] = 0;
      best_tour.count = n+1;
      *loc_best_cost_p = best_tour.cost = 
          tour_p->cost + mat[city*n + 0];
   } else if (*loc_best_cost_p > best_tour.cost) {
      *loc_best_cost_p = best_tour.cost;
   }
   pthread_mutex_unlock(&best_tour_mutex);
}  /* Check_best_tour */



/*------------------------------------------------------------------*/
void Push(tour_t* tour_p, city_t city, my_stack_t* stack_p) {
   stack_elt_t* temp = malloc(sizeof(stack_elt_t));
   temp->tour_p = Dup_tour(tour_p);
   temp->city = city;

/* pthread_mutex_lock(&(stack_p->mutex)); */

   if (stack_p->top == stack_p->allocated-1) {
      stack_elt_t** new_list = 
         malloc(2*(stack_p->allocated)*sizeof(stack_elt_t*));
      int i;
      for (i = 0; i < stack_p->allocated; i++)
         new_list[i] = stack_p->list[i];
      for (i = stack_p->allocated; i < 2*(stack_p->allocated); i++)
         new_list[i] = NULL;
      free(stack_p->list);
      stack_p->list = new_list;
      stack_p->allocated *= 2;
   }

   stack_p->top++;
   stack_p->list[stack_p->top] = temp;

/* pthread_mutex_unlock(&(stack_p->mutex)); */
}  /* Push */


/*------------------------------------------------------------------*/
tour_t* Dup_tour(tour_t* tour_p) {
   int i;
   tour_t* temp_p = malloc(sizeof(tour_t));
   temp_p->cities = malloc(n*sizeof(city_t));
   for (i = 0; i < n; i++)
      temp_p->cities[i] = tour_p->cities[i];
   temp_p->cost = tour_p->cost;
   temp_p->count = tour_p->count;
   return temp_p;
}  /* Dup_tour */

/*------------------------------------------------------------------*/
void Pop(tour_t** tour_pp, city_t* city_p, my_stack_t* stack_p) {

/* pthread_mutex_lock(&(stack_p->mutex)); */
   int top = stack_p->top;
   *tour_pp = stack_p->list[top]->tour_p;
   *city_p = stack_p->list[top]->city;
   free(stack_p->list[top]);
   stack_p->list[top] = NULL;
   stack_p->top--;
/* pthread_mutex_unlock(&(stack_p->mutex)); */
}  /* Pop */


/*------------------------------------------------------------------*/
int  Empty(my_stack_t* stack_p) {
/* pthread_mutex_lock(&(stack_p->mutex)); */
   if (stack_p->top == -1)
      return TRUE;
   else
      return FALSE;
/* pthread_mutex_unlock(&(stack_p->mutex)); */
}  /* Empty */

/*------------------------------------------------------------------*/
void Initialize_stack(my_stack_t* stack_p) {
   int i;

   stack_p->list = 
      malloc(INITIAL_STACK_SIZE*sizeof(stack_elt_t*));
   stack_p->top = -1;
   stack_p->allocated = INITIAL_STACK_SIZE;
   for (i = 0; i < stack_p->allocated; i++)
      stack_p->list[i] = NULL;
/* pthread_mutex_init(&(stack_p->mutex), NULL); */
}  /* Initialize_stack */

/*------------------------------------------------------------------*/
void Free_stack(my_stack_t* stack_p) {
   int i;

   for (i = 0; i < stack_p->allocated; i++)
      if (stack_p->list[i] != NULL)
         free(stack_p->list[i]);
   free(stack_p->list);
}  /* Free_stack */

/*------------------------------------------------------------------
 * Function:        Create_initial_records
 * Purpose:         Create the initial two-city partial tours,
 *                  and put them on the stack.
 * In arg:          my_rank
 * Out arg:         stack_p
 * Global vars in:  thread_count
 */
void Create_initial_records(long my_rank, my_stack_t* stack_p) {
   int i;
   city_t city;
   city_t cities[n+1];
   tour_t tour;
   int my_first_city, my_last_city;

   Find_initial_cities(my_rank, &my_first_city, &my_last_city);
   /* Don't call Initialize_tour since it will malloc storage for the cities */
   /* This will be lost when Push dups the tour                              */
   for (i = 0; i <= n; i++)
      cities[i] = NO_CITY;
   tour.cities = cities;
   tour.cities[0] = 0;
   tour.count = 1;
   for (city = my_first_city; city <= my_last_city; city++)  {
      tour.cost = mat[n*0 + city];
      Push(&tour, city, stack_p);
   }
}  /* Create_initial_records */

/*------------------------------------------------------------------
 * Function:        Find_initial_cities
 * Purpose:         Find the first and last cities for this thread to
 *                  use as second cities on its initial set of partial
 *                  tours
 * In args:         my_rank
 * Out args:        my_first_city_p, my_last_city_p
 * Global in vars:  n, thread_count
 */
void Find_initial_cities(long my_rank, int* my_first_city_p, 
   int* my_last_city_p) {
   int quotient, remainder;
   int my_city_count;

   quotient = (n-1)/thread_count;
   remainder = (n-1) % thread_count;
   if (my_rank < remainder) {
      my_city_count = quotient+1;
      *my_first_city_p = my_rank*my_city_count + 1;
      *my_last_city_p = *my_first_city_p + my_city_count - 1;
   } else {
      my_city_count = quotient;
      *my_first_city_p = my_rank*quotient + remainder + 1;
      *my_last_city_p = *my_first_city_p + my_city_count - 1;
   }

#  ifdef DEBUG
   printf("Thread %ld > first = %d, last = %d, count = %d\n",
      my_rank, my_first_city, my_last_city, my_city_count);
#  endif
   
}  /* Find_initial_cities */