/* Heterogeneous and Grid Computing: Assignment 2
 * MPI Implementation of the 2-dimensional ScaLAPACK
 * algorithm for matrix-multiplication
 * Author: Donal Simmie
 * Student No: 06160719
 */

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


#define MASTER 0

// function prototypes
// allocate memory for matrices
double** allocateMatrixMemory( double** matrix, int rows, int cols );
// initialise a matrix to a given value
void initialiseMatrix( double** matrix, int rows, int cols, double value );
// print matrix
void printMatrix( double** matrix, int rows, int cols );

int main( int argc, char** argv )
{
	// variabe declarations		
	int n;			// matrix dimension
	int proc_count;	// process count
	int proc_rank;	// process rank 
	int root_p;		// square root of p
	int m;			// n/r

	int i,j;		// matrix indices
	int k;			// loop index
	int ind1, ind2;	// for recalculating Bcols indices
	
	int prow, pcol;	// current processes processor grid indices
	int rank=0;

	int		**proc_grid = NULL;					// 2d processor grid
	double	**A = NULL, **B = NULL, **C = NULL;	// the matrices
	double	**Arows = NULL;	// n/r * n matrix of A data for a process row
	double  **Bcols = NULL;	// n/r * n matrix of B data for a process column (stored as horizontal blocks)

	time_t		start, end;		// for timing
	MPI_Status status;			
	MPI_Comm row_com, col_com;	// row and clumn communicators

	MPI_Init(&argc, &argv);
	
	// determine processor count and relative rank
	MPI_Comm_size (MPI_COMM_WORLD, &proc_count);
	MPI_Comm_rank (MPI_COMM_WORLD, &proc_rank );

	printf("Process %d\n", proc_rank);
	
	MPI_Finalize();

	return 0;

	root_p = sqrt(proc_count);
	n = 600;
	m = n/root_p;

	

	

	// create 2d array of processors
	proc_grid =(int **)malloc(root_p*sizeof(int *));		
	proc_grid[0] = (int *)malloc(root_p*root_p*sizeof(int));
	if(!proc_grid)
	{
		printf("Memory allocation failed \n");
		exit(1);
	}
	for(i=1; i<root_p; i++)
	{
		proc_grid[i] = proc_grid[0]+i*root_p;
		if (!proc_grid[i])
		{
			printf("Memory allocation failed \n");
			exit(1);
		}
	}
	for(i=0; i<root_p; i++)
	{
		for(j=0; j<root_p; j++)
		{
			// fill processor grid with processor ranks
			// and communicate the prow and pcol to the processor 
			// in question
			proc_grid[i][j] = rank;

			if( rank == proc_rank )
			{
				prow = i;
				pcol = j;
			}
			rank++;
			printf("%3.1d ",proc_grid[i][j]);			
		}
		printf("\n");
	}

	// create matrices
	A = allocateMatrixMemory( A, m, m );
	B = allocateMatrixMemory( B, m, m );
	C = allocateMatrixMemory( C, m, m );
	A = allocateMatrixMemory( Arows, m, m );
	B = allocateMatrixMemory( Bcols, m, m );

	// initialise matrices
	initialiseMatrix( A, m, m, 1 );
	initialiseMatrix( B, m, m, 2 );

	//start timing
	start = time(NULL);	
		
	// split communicators
	MPI_Comm_split(MPI_COMM_WORLD, prow, proc_rank, &row_com );
	MPI_Comm_split(MPI_COMM_WORLD, pcol, proc_rank, &col_com );

	// send this processes A data and receive row neighbours A data
	MPI_Allgather( A, 40000, MPI_DOUBLE, Arows, 40000, MPI_DOUBLE, row_com );
	// send this processes B data and receive column neighbours B data
	MPI_Allgather( B, 40000, MPI_DOUBLE, Bcols, 40000, MPI_DOUBLE, col_com );

	// perform local multiplication to determine this processes block of C
	for(i=0; i<m; i++)
	{
		for(j=0; j<m; j++)
		{
			C[i][j] = 0.0;
			ind1=1;
			ind2=j;
			for(k=0; k<n; k++)
			{	
				C[i][j] = C[i][j] + Arows[i][k]*Bcols[ind1][ind2];
				// recalculate indices for Bcols as they are stored horizontally
				// instead of vertically
				if( ind1 < m )
				{
					ind1++;
				}
				else
				{
					ind1=1;
					ind2+=m;
				}
			}
		}
	}// end multiplication

	// wait for everyone to finish
	//MPI_Barrier( row_com );
	//MPI_Barrier( row_com );

	end = time(NULL);//end timing
		
	if( proc_rank == MASTER )
	{
		// print time taken - which has some communication overhead
		// but this is negligible when dealing with large n
		printf("\nC[m-1][m-1]: %1.1f\n",C[m-1][m-1]);
		//printf("\nC(%d)\n",proc_rank);
		//printMatrix(C,m,m);
		printf("Elapsed time = %3d seconds.\n", end - start);			
	}

	// MPI_Finalize is the last instruction of the program
	MPI_Finalize();

	// release allocated memory
	free(A[0]);
	free(A);
	free(B[0]);
	free(B);
	free(C[0]);
	free(C);
	free(Arows[0]);
	free(A);
	free(Bcols[0]);
	free(B);

	return 0;
}

double** allocateMatrixMemory( double** matrix, int rows, int cols )
{	
	int i;
	matrix =(double **)malloc(rows*sizeof(double *));		 // rows
	matrix[0] = (double *)malloc(rows*cols*sizeof(double));  // data
	if(!matrix)
	{
		printf("Memory allocation failed \n");
		exit(1);
    }
	for(i=1; i<rows; i++)
    {
		matrix[i] = matrix[0]+i*cols;
		if (!matrix[i])
		{
			printf("Memory allocation failed \n");
			exit(1);
		}
    }
	return matrix;
}

void initialiseMatrix( double** matrix, int rows, int cols, double value )
{
	int i,j;
	
	for(i=0;i<rows;i++)
	{
		for(j=0;j<cols;j++)
		{
			matrix[i][j] = value;
		}
	}
}

void printMatrix( double** matrix, int rows, int cols )
{
	int i,j;

	for(i=0;i<rows;i++)
	{
		for(j=0;j<cols;j++)
		{
			printf("%1.0f ", matrix[i][j] );
		}
		printf("\n");
	}
	printf("\n");
}
// end file Assignment2.c