/******************************************************************************
 * Copyright (C) International Computer Science Institute, 1995.  COPYRIGHT   *
 * NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject  *
 * to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in     *
 * the file "Doc/License" of the Sather distribution.  The license is also    *
 * available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA.   *
 * -----> Please email comments to "sather-bugs@icsi.berkeley.edu". <-------- *
 *                                                                            *
 * pmatrix.c                                                                  *
 * implements naive parallel matrix multiplication used to 		      *
 * compare pSather and C       						      *
 * Version 1.0, June 1996, by Claudio Fleiner		                      *
 ******************************************************************************/

#include "am.h"
#include "times.h"

long next() 
{
	static long seed=997;
	seed=seed*123457+97;
	if(seed<0) seed= -seed;
	return seed;
}

float rnd_floating() { return (float)(next()%1000000)/1000000.0; }

struct MATRIX {
	int rows;
	int cols;
	float *t;
};
#define M(m,x,y) m.t[x+y*m.cols]

struct MATRIX new_matrix(int cols,int rows)
{
	struct MATRIX t;
	t.rows=rows;
	t.cols=cols;
	t.t=(float *)malloc(sizeof(float)*cols*rows);
	return t;
}

void print_matrix(struct MATRIX t) {
	int i,j;
	for(i=0;i<t.rows;i++) {
		printf("| ");
		for(j=0;j<t.cols;j++) 
			printf("%f ",M(t,j,i));
		printf("|\n");
	}
	printf("\n");
}

struct MATRIX a,b,r;
sema_t sm;
void matrix_mult_t(vnn_t from,int l,int u) {
	int x,y,z;
	int r_cols,a_cols,b_cols,a_rows;
	float *rt,*at,*bt;

	printf("thread working on %d - %d\n",l,u);
	for(x=l;x<u;x++)
		for(y=0;y<b.cols;y++) {
			M(r,y,x)=0.0;
			for(z=0;z<a.cols;z++) {
				M(r,y,x)+=M(a,z,x)*M(b,y,z);
			}
		}
	SEMA_SIGNAL(sm);
}


int nthr=4;
struct MATRIX matrix_mult(struct MATRIX a,struct MATRIX b) {
	int i;
	sm=SEMA_CREATE(0);
	r=new_matrix(a.rows,a.rows);
	for(i=0;i<nthr;i++) thr_create_2(AM_MY_PROC,(handler_2_t)matrix_mult_t,i*a.rows/nthr,(i+1)*a.rows/nthr);
	for(i=0;i<nthr;i++) SEMA_WAIT(sm);
	return r;
}


int main(int argc,char *argv[]) {
	int i,j;
	int size=atoi(argv[1]);
	struct MATRIX c;
	TIMES t;
	/*
	a=new_matrix(1,2);
	b=new_matrix(2,1);
	M(a,0,0)=1;
	M(a,0,1)=2;
	M(b,0,0)=3;
	M(b,1,0)=4;
	print_matrix(a);
	print_matrix(b);
	c=matrix_mult(a,b);
	print_matrix(c);
	c=matrix_mult(b,a);
	print_matrix(c);
	*/

	am_enable(1,argc,argv);
	a=new_matrix(size,size);
	b=new_matrix(size,size);
	thr_setconcurrency(atoi(argv[2]));
	nthr=atoi(argv[3]);
	for(i=0;i<size;i++)
		for(j=0;j<size;j++) {
			M(a,i,j)=rnd_floating();
			M(b,i,j)=rnd_floating();
		}

	t=times_now();
	c=matrix_mult(a,b);
	t=times_elapsed(t);
	am_disable();
	printf("%s",times_str(t));
/*	print_matrix(c); */
}