#include "stdio.h"
#include "mpi.h"	/* This brings in pre-defined MPI constants, ... */

void asemble(int at[][], int ml, int nl, MPI_Comm comm, int b[][], int m, int n, int p);

void main(int argc, char *argv[]) 
{
	int n, m, nv, nl, mv, ml, i, il, iv, j, jl, jv;
	int p, ndim, reorder, ierr;
	int master, me, Iam, source, dest, tag;
	int dims[2], coord[2];
	int period[2];
	int a[3][2], at[2][3], b[4][9];
	MPI_Status status;
	MPI_Request req;
	MPI_Comm grid_comm;

	/* Starts MPI processes ... */
	MPI_Init(&argc, &argv);					/* starts MPI */
	MPI_Comm_rank(MPI_COMM_WORLD, &Iam);	/* get current process id */
	MPI_Comm_size(MPI_COMM_WORLD, &p);		/* get number of processes */

	master = 0;	/* 0 is defined as the master processor */
	period[0] = 0; period[1] = 0;	/* no cyclic boundary in either index */
	tag = 0;	/* a tag is not required in this case, set it to zero */
	dest = 0;	/* results are sent back to master */

	n = 9; m = 4; nv = 3; mv = 2;
	nl = n/nv; ml = m/mv;
	ndim = 2; reorder = 1;

	/* create cartesian topology for matrix */
	dims[0] = nv;
	dims[1] = mv;
	MPI_Cart_create(MPI_COMM_WORLD, ndim, dims, period, reorder, &grid_comm);
	MPI_Comm_rank(grid_comm, &me);
	MPI_Cart_coords(grid_comm, me, ndim, coord);
	iv = coord[0];
	jv = coord[1];

	/* define local matrix according to virtual grid coordinates, (iv,jv) */
	for (il=0; il<nl; il++)
	{
		for (jl=0; jl<ml; jl++)
		{
			i = il + iv*nl;
			j = jl + jv*ml;
			a[il][jl] = i*10 + j;
		}
	}

	printf("%d: Before Transpose:\n", Iam);
	for (i=0; i<nl; i++)
	{
		for (j=0; j<ml; j++)
		{
			printf("%5d", a[i][j]);
		}
		printf("\n");
	}

	/* perform transpose on local matrix */
	for (il=0; il<nl; il++)
	{
		for (jl=0; jl<ml; jl++)
		{
			at[jl][il] = a[il][jl];
		}
	}

	/* send "at" to Master for asembly and printing */
	MPI_Isend(at, ml*nl, MPI_INT, master, tag, grid_comm, &req);

	/* Master asembles all local transposes into final matrix and print */
	if(Iam == master) {
		asemble(at, ml, nl, grid_comm, b, m, n, p);
		MPI_Wait(&req, &status);	/* make sure all sends done */
	}
	MPI_Finalize();		/* let MPI finish up ... */
}

void asemble(int at[2][3], int ml, int nl, MPI_Comm comm, int b[4][9], int m, int n, int p)
{
	int tag, source, ierr, ndim;
	int iv, jv, i, j, il, jl, coord[2];
	MPI_Status status;
	tag = 0;

	/* The Master asembles the final (transposed) matrix from local copies and print */
	for (source=0; source<p; source++)
	{
		MPI_Cart_coords(comm, source, ndim, coord);
		MPI_Recv(at, ml*nl, MPI_INT, source, tag, comm, &status);
		iv = coord[0];
		jv = coord[1];
		for (jl=0; jl<nl; jl++)
		{
			j = jl + iv*nl;		/* swap iv and jv for transpose */
			for (il=0; il<ml; il++)
			{
				i = il + jv*ml;
				b[i][j] = at[il][jl];
			}
		}
	}

	printf("\nAfter Transpose:\n");
	for (i=0; i<m; i++)
	{
		for (j=0; j<n; j++)
		{
			printf("%5d", b[i][j]);
		}
		printf("\n");
	}
}