/* File:     
 *     ex4.2_pth_mat_vect.c
 *
 * Purpose:  
 *     Computes a parallel matrix-vector product. The matrix
 *     is physically distributed by block rows.  The vector x is
 *     logically distributed by blocks.  The vector y is 
 *     physically distributed by blocks. This is a modified
 *     version of pth_mat_vec.c.
 *
 * Input:
 *     m, n: order of matrix
 *     x, A: the vector and the matrix to be multiplied
 *
 * Output:
 *     y: the product vector
 *
 * Compile:  
 *     gcc -g -Wall -o ex4.2_pth_mat_vect ex4.2_pth_mat_vect.c -lpthread
 * Usage:
 *     pth_mat_vect <thread_count>
 *
 * Notes:  
 *     1.  Storage for A, x, y is dynamically allocated.
 *     2.  Number of threads (thread_count) should evenly divide both 
 *         m and n.  The program doesn't check for this.
 *     3.  We use a 1-dimensional array for A and compute subscripts
 *         using the formula A[i][j] = A[i*n + j]
 *     4.  The vector x is globally shared.  A and y are distributed
 *         among the the threads.
 *
 * IPP:    Exercise 4.2 (p. 200)
 */

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <semaphore.h>
#include "timer.h"

/* Global variables */
int     thread_count;
int     m, n;
double* x;

/* Semaphores */
sem_t* sem_a;
sem_t* sem_y;

/* Mutexs */
pthread_mutex_t bar_mutex, max_mutex;
pthread_cond_t bar_cond;
int bar_count = 0;
double elapsed = 0;

/* Serial functions */
void Usage(char* prog_name);
void Read_matrix(char* prompt, double A[], int m, int n);
void Read_vector(char* prompt, double x[], int n);
void Print_matrix(char* title, double A[], int m, int n);
void Print_vector(char* title, double y[], double m);
void Semaphores_wait(sem_t* sem, long rank);
void Semaphores_post(sem_t* sem, long rank);
void My_barrier(void);

/* Parallel function */
void *Pth_mat_vect(void* rank);

/*------------------------------------------------------------------*/
int main(int argc, char* argv[]) {
    long       thread;
    pthread_t* thread_handles;
    int tmp1, tmp2;
    
    if (argc != 2) Usage(argv[0]);
    thread_count = strtol(argv[1], NULL, 10);
    /* Allocate arrays of semaphores to schedule input */
    sem_a = malloc(thread_count*sizeof(sem_t));
    sem_y = malloc(thread_count*sizeof(sem_t));
    
    /* Initialize semaphores */
    for(thread = 0; thread < thread_count; thread++) {
        if(thread == 0) {
            tmp1 = sem_init(&sem_a[thread], 0, 1);
            tmp2 = sem_init(&sem_y[thread], 0, 0);
        }
        else {
            tmp1 = sem_init(&sem_a[thread], 0, 0);
            tmp2 = sem_init(&sem_y[thread], 0, 0);
        }
        if(tmp1 == -1 || tmp2 == -1) {
            fprintf(stderr, "sem_init failed for thread %ld\n", thread);
            free(sem_a);
            free(sem_y);
            exit(-1);
        }
    }

    /* Allocate array for threads handles*/
    thread_handles = malloc(thread_count*sizeof(pthread_t));
    
    /* Initialize the mutexes */
    pthread_mutex_init(&bar_mutex, NULL);
    pthread_mutex_init(&max_mutex, NULL);
    pthread_cond_init(&bar_cond, NULL);
    
    printf("Enter m and n\n");
    scanf("%d%d", &m, &n);

    x = malloc(n*sizeof(double));

    Read_vector("Enter the vector x", x, n);
    
    printf("Enter the matrix\n");
    
    /* Start threads */
    for (thread = 0; thread < thread_count; thread++)
        pthread_create(&thread_handles[thread], NULL,
            Pth_mat_vect, (void*) thread);
    /* Wait for threads to complete */
    for (thread = 0; thread < thread_count; thread++)
        pthread_join(thread_handles[thread], NULL);

    printf("Elapsed time = %e\n", elapsed);

    for(thread = 0; thread < thread_count; thread++) {
        sem_destroy(&sem_a[thread]);
        sem_destroy(&sem_y[thread]);
    }

    pthread_mutex_destroy(&bar_mutex);
    pthread_mutex_destroy(&max_mutex);
    pthread_cond_destroy(&bar_cond);
    
    free(thread_handles);
    free(sem_a);
    free(sem_y);
    free(x);

   return 0;
}  /* main */


/*------------------------------------------------------------------
 * Function:  Usage
 * Purpose:   print a message showing what the command line should
 *            be, and terminate
 * In arg :   prog_name
 */
void Usage (char* prog_name) {
   fprintf(stderr, "usage: %s <thread_count>\n", prog_name);
   exit(0);
}  /* Usage */

/*------------------------------------------------------------------
 * Function:    Read_matrix
 * Purpose:     Read in the matrix
 * In args:     prompt, m, n
 * Out arg:     A
 */
void Read_matrix(char* prompt, double A[], int m, int n) {
   int             i, j;

   printf("%s\n", prompt);
   for (i = 0; i < m; i++) 
      for (j = 0; j < n; j++)
         scanf("%lf", &A[i*n+j]);
}  /* Read_matrix */

/*------------------------------------------------------------------
 * Function:        Read_vector
 * Purpose:         Read in the vector x
 * In arg:          prompt, n
 * Out arg:         x
 */
void Read_vector(char* prompt, double x[], int n) {
   int   i;

   printf("%s\n", prompt);
   for (i = 0; i < n; i++) 
      scanf("%lf", &x[i]);
}  /* Read_vector */

/*------------------------------------------------------------------
 * Function:        Semaphores_wait
 * Purpose:         wait to get a lock
 * In arg:          sem, rank
 */
void Semaphores_wait(sem_t* sem, long rank) {
    if(sem_wait(sem) == -1) {
        fprintf(stderr, "Sem_wait failed for rank %ld\n", rank);
        exit(0);
    }
}  /* Semaphores_wait */

/*------------------------------------------------------------------
 * Function:        Semaphores_post
 * Purpose:         release semaphore lock
 * In arg:          sem, rank
 */
void Semaphores_post(sem_t* sem, long rank) {
    if(sem_post(sem) == -1) {
        fprintf(stderr, "Sem_post failed for rank %ld\n", rank);
        exit(0);
    }
}  /* Semaphores_post */


/*------------------------------------------------------------------
 * Function:       Pth_mat_vect
 * Purpose:        Multiply an mxn matrix by an nx1 column vector
 * In arg:         rank
 * Global in vars: A, x, m, n, thread_count
 * Global out var: y
 */
void *Pth_mat_vect(void* rank) {
    long my_rank = (long) rank;
    int local_i, j;
    int local_m = m/thread_count; 
    
    double* local_A = malloc(local_m*n*sizeof(double));
    double* local_y = malloc(local_m*sizeof(double));
    
    double start, finish, my_elapsed;
    
    /* Serialize input of local_A using semaphores */
    Semaphores_wait(&sem_a[my_rank], my_rank);

    for (local_i = 0; local_i < local_m; local_i++) 
        for (j = 0; j < n; j++)
            scanf("%lf", &local_A[local_i*n+j]);
    
    if(my_rank+1 < thread_count)
        Semaphores_post(&sem_a[my_rank+1], my_rank+1);
    
    /* record timer */
    My_barrier();
    GET_TIME(start);
    for (local_i = 0; local_i < local_m; local_i++) {
        local_y[local_i] = 0.0;
        for (j = 0; j < n; j++)
            local_y[local_i] += local_A[local_i*n + j]*x[j];
    }
    GET_TIME(finish);
    my_elapsed = finish-start;

    pthread_mutex_lock(&max_mutex);
    if (my_elapsed > elapsed) elapsed = my_elapsed;
    pthread_mutex_unlock(&max_mutex);
        
    /* Serialize output of local_y using semaphores */
    if (my_rank + 1 == thread_count)
        Semaphores_post(&sem_y[0], 0);
    
    Semaphores_wait(&sem_y[my_rank], my_rank);
    
    if(my_rank == 0) 
        Print_vector("The product is\n", local_y, local_m);
    else 
        Print_vector("", local_y, local_m);
    fflush(stdout);
    
    if (my_rank+1 < thread_count)
        Semaphores_post(&sem_y[my_rank+1], my_rank+1);
    else 
        printf("\n");

    free(local_A);
    free(local_y);
    
    return NULL;
}  /* Pth_mat_vect */

/*--------------------------------------------------------------*/

void My_barrier(void) {
    pthread_mutex_lock(&bar_mutex);
    bar_count++;
    if (bar_count == thread_count) {
        bar_count = 0;
        pthread_cond_broadcast(&bar_cond);
    } else {
        while (pthread_cond_wait(&bar_cond, &bar_mutex) != 0);
    }
    pthread_mutex_unlock(&bar_mutex);
}  /* My_barrier */

/*------------------------------------------------------------------
 * Function:    Print_matrix
 * Purpose:     Print the matrix
 * In args:     title, A, m, n
 */
void Print_matrix( char* title, double A[], int m, int n) {
   int   i, j;

   printf("%s\n", title);
   for (i = 0; i < m; i++) {
      for (j = 0; j < n; j++)
         printf("%4.1f ", A[i*n + j]);
      printf("\n");
   }
}  /* Print_matrix */


/*------------------------------------------------------------------
 * Function:    Print_vector
 * Purpose:     Print a vector
 * In args:     title, y, m
 */
void Print_vector(char* title, double y[], double m) {
   int   i;

   printf("%s", title);
   for (i = 0; i < m; i++)
      printf("%4.1f ", y[i]);
   //printf("\n");
}  /* Print_vector */