/**********************************************************************
		 C-DAC Tech Workshop : hyPACK-2013
                             October 15-18,2013


 Example 1.3      : mpi-Pthreads-infinity-norm.c

 Objective        : Write an MPI-Pthreads program to calculate Infinity 
                    norm of a matrix using row wise Block-striped 
                    partitioning.  

                    Pthreads and MPI Library calls are used.

                    This Example demonstrates the use of:

                    pthread_create()
                    pthread_join()
                    pthread_mutex_lock()
                    pthread_mutex_unlock()

                    MPI_Init
                    MPI_Comm_rank
                    MPI_Comm_size
                    MPI_Barrier
                    MPI_Bcast
                    MPI_Allgather
                    MPI_Reduce
                    MPI_Finalize

 Input             : The input file holding the matrix.

 Output            : Infinity Norm of the given Matrix.

 Created           :MAY-2012 
                         
 E-mail      	   : betatest@cdac.in                                          


************************************************************************/

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

int  MyRank, currentRow, MyNoofRows, NoofCols, GlobalIndex = -1;
int  flag = 0, maxflag = 0, rowlimit;
float **Matrix, *Vector, max;

pthread_mutex_t mutex_Row = PTHREAD_MUTEX_INITIALIZER;
pthread_mutex_t mutex_Flag = PTHREAD_MUTEX_INITIALIZER;
pthread_mutex_t mutex_Max = PTHREAD_MUTEX_INITIALIZER;

/* Routine executed by each thread */

void * MyPartOfCalc(int Id)
{

	int    myRow, col;
	float  sum;

	if (flag == 0) 
        {
		pthread_mutex_lock(&mutex_Flag);
		rowlimit = currentRow + MyNoofRows;
		flag++;
		pthread_mutex_unlock(&mutex_Flag);
	}

	while (1) 
        {

		/*
		 * Thread selects the row of Matrix on which it has to do the
		 * operation
		 */

		pthread_mutex_lock(&mutex_Row); 
                {
			if (currentRow >= rowlimit) 
                        {
				pthread_mutex_unlock(&mutex_Row);
				pthread_exit(0);
			}
			myRow = currentRow;
			currentRow++;
		}
		pthread_mutex_unlock(&mutex_Row);

		/*
		 * Perform the addition  on the row selected and and store in
		 * max if it is the maximum value until now out of computed
		 * sums
		 */

		printf(" Thread Id %d of process with Rank %d operated on Matrix Row %d\n", Id, MyRank, myRow);
		sum = 0.0;

		for (col = 0; col < NoofCols; col++)
		sum += fabs(Matrix[myRow][col]);

		pthread_mutex_lock(&mutex_Max);

		if (maxflag == 0) 
                 {
			max = sum;
			maxflag = 1;
		  } 
                else if (sum > max)
			max = sum;

		pthread_mutex_unlock(&mutex_Max);
	}

}


main(int argc, char **argv)
{

	int   iproc, irow, icol, modval, divval, *Displacement;
	int  iprocb, *ArrayNoofRows, Numprocs, Root = 0, NoofRows, VectorSize;
	float  Result;

	pthread_t  *threads;
        FILE  *fp;


	/* MPI Initialisation ... */

	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &MyRank);
	MPI_Comm_size(MPI_COMM_WORLD, &Numprocs);

	/* Validity checking for minimum number of processors */

	if (Numprocs < 2) 
        {
		printf("Invalid Number of Processors ..... \n");
		printf("Numprocs must be greater than 1 ......\n");
		MPI_Finalize();
		exit(0);
	 }

        /* Read the Matrix from the file. */


        fp = fopen("Norm_Data.dat", "r");
        if (!fp)
        {
          printf("\nUnable to open the file Norm_Data.dat");
          exit(0);
        }

        fscanf(fp, "%d", &NoofRows);
        fscanf(fp, "%d", &NoofCols);

        printf("\n Rows: %d  Col:%d.", NoofRows, NoofCols);

	MPI_Barrier(MPI_COMM_WORLD);

	MPI_Bcast(&NoofRows, 1, MPI_INT, Root, MPI_COMM_WORLD);
	MPI_Bcast(&NoofCols, 1, MPI_INT, Root, MPI_COMM_WORLD);

	/* Validity checking for negative sizes of Matrix */

	if (NoofRows < 1 || NoofCols < 1) 
        {
		printf("The number of rows or columns or size of Vector should be atleast one\n");
		MPI_Finalize();
		exit(-1);
	}
	/* Validity checking for minimum number of Rows of Matrix */

	if (NoofRows < Numprocs) 
        {
		printf("The number of rows of Matrix should be greater than number of processors\n");
		MPI_Finalize();
		exit(-1);
	}
	/* Allocating and Populating the Matrix */

	Matrix = (float **) malloc(NoofRows * sizeof(float *));
	for (irow = 0; irow < NoofRows; irow++)
	Matrix[irow] = (float *) malloc(NoofCols * sizeof(float));

	for (irow = 0; irow < NoofRows; irow++)
	for (icol = 0; icol < NoofCols; icol++)
        fscanf(fp, "%f", &Matrix[irow][icol]);


	/* Storing the number of Rows to be operated by each process in array */

	modval = NoofRows % Numprocs;
	divval = NoofRows / Numprocs;
	MyNoofRows = (MyRank < modval ? divval + 1 : divval);

	ArrayNoofRows = (int *) malloc(Numprocs * sizeof(int));
	MPI_Allgather(&MyNoofRows, 1, MPI_INT, ArrayNoofRows, 1, MPI_INT, MPI_COMM_WORLD);

	/* Storing the starting Row to be  operated by each process in array */

      Displacement = (int *) malloc(Numprocs * sizeof(int));
      Displacement[0] = 0;
      for (iproc = 1; iproc < Numprocs; iproc++)
      Displacement[iproc] = Displacement[iproc - 1] + ArrayNoofRows[iproc - 1];

	currentRow = Displacement[MyRank];

	MPI_Barrier(MPI_COMM_WORLD);

	/*
	 * Call threads equal to number of Rows to be processed by this
	 * process
	 */

	threads = (pthread_t *) malloc(sizeof(pthread_t) * MyNoofRows);

	for (irow = 0; irow < MyNoofRows; irow++)
	pthread_create(&threads[irow], NULL, (void *(*) (void *)) MyPartOfCalc, (void *) irow);

	MPI_Barrier(MPI_COMM_WORLD);
	for (irow = 0; irow < MyNoofRows; irow++)
	pthread_join(threads[irow], NULL);

	MPI_Barrier(MPI_COMM_WORLD);

	MPI_Reduce(&max, &Result, 1, MPI_FLOAT, MPI_MAX, Root, MPI_COMM_WORLD);


	if (MyRank == Root) 
        {
		printf("The matrix is :\n\n");
		for (irow = 0; irow < NoofRows; irow++) 
                  {
			for (icol = 0; icol < NoofCols; icol++)
			printf("%.1f ", Matrix[irow][icol]);
			printf("\n");
		  }
		printf("The infinity norm is %f \n", Result);
	}
	MPI_Finalize();

}