/* 
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-26
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

 A simple pthreads parallel matrix multiplication program
 (Matrix_A  X  Matrix_B) => Matrix_C

*/ 

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>

/* Global storage for matrices A,B,C */
int **mA,**mB,**mC;

/* Storage for anything I need to pass to the worker threads */ 
typedef struct {
  int thread_number; 
  int size;                       /* Total size of each matrix */  
  int row_start;  
  int row_end; 
} tpackage_t; 


/* Worker routine for each thread: 

   This rouine will compute the rows of matrices C
   from [row_start, row_end).
*/ 
void *mult_worker( void *arg ){
  
  int thread_number, row_start, row_end, size; 
  int i,j,k,status; 
  tpackage_t *p = (tpackage_t*) arg; 

  status = pthread_detach( pthread_self() ); 
  if( status != 0 ){
    fprintf(stderr,"status != 0");
    exit(2);
  }

  /* Extract this threads "package" */ 
  thread_number = p->thread_number; 
  row_start = p->row_start; 
  row_end = p->row_end; 
  size = p->size; 

  printf( "I AM thread # %4d ... let me work!\n", thread_number );
   
  /* Initialize these rows of mC: */ 
  for( i=row_start; i < row_end; i++ ){
    for( j=0; j < size; j++ ){
      mC[i][j] = 0; 
    }
  }

  /* Compute the product matric mC: */ 
  for( i=row_start; i < row_end; i++ ){
    for( j=0; j < size; j++ ){
      for( k=0; k < size; k++ ){
	mC[i][j] += mA[i][k]*mB[k][j];
      }
    }
  }

  free(p); 

  printf( "I AM thread # %4d ... just finished!\n", thread_number );
}


int 
main( int argc, char* argv[]){

  int n, P, i, j, ti, num_strip, row_start, row_end; 
  pthread_t *threads;
  tpackage_t *p;                  /* Argument list each thread will get */ 
  
  /* Argument checking */
  if( argc != 3 ){
    fprintf( stderr, "Expecting three input arguments\n" ); 
    exit(1); 
  }
  /* Input the size of each matrix: */ 
  n = atoi(argv[1]);
  /* Input the number of processes/theads */
  P = atoi(argv[2]); 
  
  /* Create matrix A */ 
  mA = (int**) malloc( n*sizeof(int*) ); 
  if( !mA ){
    fprintf(stderr,"Memory couldn't be created\n"); 
    exit(1);
  }
  for( i=0; i < n; i++ ){
    mA[i] = malloc( n*sizeof(int) ); 
    if( !mA[i] ){ 
      fprintf(stderr,"Memory couldn't be created\n"); 
      exit(1);
    }
  }

  /* Create matrix B */ 
  mB = (int**) malloc( n*sizeof(int*) ); 
  if( !mB ){
    fprintf(stderr,"Memory couldn't be created\n"); 
    exit(1);
  }
  for( i=0; i < n; i++ ){
    mB[i] = malloc( n*sizeof(int) ); 
    if( !mB[i] ){ 
      fprintf(stderr,"Memory couldn't be created\n"); 
      exit(1);
    }
  }

  /* Create matrix C */ 
  mC = (int**) malloc( n*sizeof(int*) ); 
  if( !mC ){
    fprintf(stderr,"Memory couldn't be created\n"); 
    exit(1);
  }
  for( i=0; i < n; i++ ){
    mC[i] = malloc( n*sizeof(int) ); 
    if( !mC[i] ){ 
      fprintf(stderr,"Memory couldn't be created\n"); 
      exit(1);
    }
  }
  /* Create a number of pthread_t objects: */ 
  threads = (pthread_t *) malloc( P*sizeof(pthread_t) ); 

  /* Initialize matrix A: */ 
  for( i=0; i < n; i++ ){
    for( j=0; j < n; j++ ){
      mA[i][j] = 1; 
    }
  }
  /* Initialize matrix B: */ 
  for( i=0; i < n; i++ ){
    for( j=0; j < n; j++ ){
      mB[i][j] = 1; 
    }
  }

  /* Divide work up by rows and crate threads to do the work */ 
  num_strip = (int) (n/P);
  for( ti=0; ti < P; ti++ ){
    row_start = 0 + num_strip*ti; 
    /* The last thread will do any extra rows if needed: */ 
    if( ti < P-1 ){
      row_end = row_start + num_strip; 
    }else{
      row_end = n;
    }
    printf( "Thread %3d will work on rows [%6d,%6d), doing a total of %6d rows\n", 
	    ti, row_start, row_end, row_end-row_start ); 
    
    /* Create and fill the args needed for the worker threads: */ 
    p = (tpackage_t*) malloc( sizeof( tpackage_t ) ); 
    p->thread_number = ti; 
    p->row_start = row_start; 
    p->row_end = row_end; 
    p->size = n; 

    pthread_create( &threads[ti], NULL, mult_worker, (void *) p ); 
  }

  /* Synchronize the completion of all threads */
  for( ti=0; ti < P; ti++ ){
    pthread_join(threads[ti],NULL); 
    printf( "Thread %d is joining the fray\n", ti ); 
  }
  /* Clear the thread memory ... no longer needed */ 
  free(threads);

  /* Display the resulting product matrix: */ 
  for(i = 0; i < n; i ++) {
    for (j = 0; j < n; j++) {
      printf("%5d ",mC[i][j]);
    }
    printf("\n");
  }


}