/* Appendix A. Source Code Listing. 					*/
/*************************************************************************
**	Project : CIS 629 project #1 Cache Memories.			**
**	Author 	: Sameer Shende 					**
**	Location: /home/grads/sameer/cis629/p10				**
**	File	: multiply.m						**
**	Description : This program runs on MASPAR-1 and calculates the  **
**		      product of two matrices C = A x B 		**
*************************************************************************/

#include <stdio.h>
#include <mpl.h>


int  init_matrix_a(); /* initialise the a matrix */
int  init_matrix_b(); /* initialise the b matrix */
void get_matrix_c();  /* initialise the c matrix */
void show_matrix_c(); /* display the C matrix 	 */
void check_matrix();  /* verify the matrix */
void rearrange_a_and_b(); /* part of multiply */
void good_method();   /* for verification. FROM MASPAR MANUAL */
void matrix_multiply(); /* THE MATRIX MULTIPLICATION ROUTINE FOR SUBMISSION */

/* Timer routines */
void dpuTimerStart();
unsigned long dpuTimerTicks();
double dpuTimerConst();

plural int c, a, b; /* Matrices a b and c. c = a x b */
int A[64][64]; /* A on FE */
int B[64][64]; /* B on FE */
int C[64][64]; /* C = A x B */

main()
{
	init_matrix_a();	/* initialise a */	
	init_matrix_b();	/* initialise b */	

	dpuTimerStart(); 	/* Timer call */
	matrix_multiply(); 	/* Multiply call */
  	printf("MATRIX MULTIPLICATION 64x64 DPU cycles: %d\n", 
		dpuTimerTicks());
  	printf("MATRIX MULTIPLICATION 64x64 DPU time:   %f\n", 
		(float)dpuTimerTicks() * (float)dpuTimerConst());


	/*
	good_method(); */ /* For verification */
	get_matrix_c(); /* get C on front end. Transfers c-> C */
/*
	show_matrix_c();
*/

}

void get_matrix_c(void)
{ /* transfers c from DPU to FE */
	int i,j,no;
	/* on the front end */

	for(i=0;i<nxproc;i++)
		for(j = 0;j< nyproc; j++)
		{
			no = i*nxproc + j;
			C[i][j] = proc[no].c;
		}
	
}

void get_matrix_a(void)
{ /* transfers a from DPU to FE */
	int i,j,no;
	/* on the front end */

	for(i=0;i<nxproc;i++)
		for(j = 0;j< nyproc; j++)
		{
			no = i*nxproc + j;
			A[i][j] = proc[no].a;
		}
	
}

void get_matrix_b(void)
{ /* transfers b from DPU to FE */
	int i,j,no;
	/* on the front end */

	for(i=0;i<nxproc;i++)
		for(j = 0;j< nyproc; j++)
		{
			no = i*nxproc + j;
			B[i][j] = proc[no].b;
		}
	
}

int init_matrix_a(void)
{ /* Put entries in a */
	a = ixproc* 3 + iyproc;
}

int init_matrix_b(void)
{ /* Put entries in b */
	b = ixproc + iyproc;
}

void rearrange_a_and_b() 
{ /* Arrange a and b so that the matrix can be fed in */
	int i;
		
	for(i=0;i<nxproc+1;i++) /* does n+1 times */
	{
		if (i >= iyproc)
			xnetE[1].a = a;
		if (i >= ixproc)
			xnetS[1].b = b;
	}

}

void show_matrix_c()
{ /* displays C */
int i,j;

	for(i=0;i<nxproc;i++)
		for(j = 0;j< nyproc; j++)
		{
			printf("C[%d][%d] = %d\n", i,j,C[i][j]);
		}
	
}

void show_matrix_a()
{ /* displays A */
int i,j;

	for(i=0;i<nxproc;i++)
		for(j = 0;j< nyproc; j++)
		{
			if(((i < 16)&&(j < 16))|| (i== nxproc - 1))
			printf("A[%d][%d] = %d\n", i,j,A[i][j]);
		}
	
}

void show_matrix_b()
{ /* displays b */
int i,j;

	for(i=0;i<nxproc;i++)
		for(j = 0;j< nyproc; j++)
		{
			if((i < 16)&&(j < 16))
			printf("B[%d][%d] = %d\n", i,j,B[i][j]);
		}
	
}

/* This piece of code has been taken from the MASPAR programming manual
verbatim - but it was used to check the results of matrix_multiply only.
The results match for all cases. This is more efficient. Reference : This 
method appears in MASPAR Programming Language (Ansi C Compatible MPL)
User Guide software version 3.1 in Section 3-10 (Data Parallel Algorithms)
under More Examples. This is in /usr/maspar/doc/MPL_ug.ps */

void good_method() /* from MASPAR manual */
{
	int i;
	for (i=0;i<nxproc-1;i++)
	{
		if (iyproc > i)
			xnetW[1].a = a;
		if (ixproc > i)
			xnetN[1].b = b;
	}

	c = 0;
	for (i=0;i<nxproc;i++)	
	{
		c += a*b;
		xnetW[1].a = a;
		xnetN[1].b = b;
	}
}
	
			
/* THE MAIN MATRIX MULTIPLICATION ROUTINE - THIS IS THE SUBMISSION */
/* To implement this concepts were taken from systolic arrays and the
blocked mesh matrix multiplication method. Reference : "The design and
Analysis of Parallel Algorithms" - Selim G Akl. Prentice Hall 1989.
page 180 section 7.3.1 "Mesh multiplication". This section describes how
to implement the matrix multiplication on a mesh by staggering A and B.
This implementation takes 3n-2 steps as it waits for the first n elements
to arrive, then P(n,n) has to wait n-1 steps for the first element to arrive
and then has to wait for n-1 steps for the last element to finish. This is
not the most optimum style and can be improved further. */

void matrix_multiply() /* from the book */
{
	int iteration =0;
	plural int count = 0;

	rearrange_a_and_b(); /* stagger a and b as required */

	c = 0;
	for(iteration=0;iteration < 3*nxproc - 2; iteration ++)
	{
		if ((iteration >= ixproc)&&(iteration >= iyproc))
		{
			count ++;
			if (count <= nxproc)
				c += a*b;
/*
			if((iproc < 5)||(iproc == 63*64+1) ||(iproc == 4095))
			{
				printf("it = %d************\n", iteration);
				p_printf("iproc = %d, a = %d, b = %d, c=%d\n", iproc, a, b, c);
			}
*/
		}
		/* shift the elements */
		xnetE[1].a = a;
		xnetS[1].b = b;
	}

}

/* EOF multiply.m */
