/******************************************************************************
 * 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". <-------- *
 *                                                                            *
 * matrix.cc                                                                  *
 * implements naive matrix multiplication used to compare Sather and C++      *
 * Version 1.0, June 1996, by Claudio Fleiner		                      *
 ******************************************************************************/
#include <stream.h>
#include "rnd.h"
#include "times.h"

template <class T>
class MATRIX {
	T* t;
    public:
	int rows,cols;
	MATRIX() {
		cols=0;
		rows=0;
		t=0;
	}
	MATRIX(int c,int r) {
		rows=r;
		cols=c;
		t=new T[r*c];
	}

	inline T &operator()(int x,int y) {
		return t[x+y*cols];
	}

	MATRIX<T>::operator*(MATRIX<T> &a) {
		MATRIX<T> r(rows,rows);
		int a_cols=a.cols;
		int r_cols=r.cols;
		int self_cols=cols;
		int solf_rows=rows;
		float *rt=r.t;
		float *at=a.t;
		float *selft=t;
		for(int x=0;x<rows;x++)
			for(int y=0;y<rows;y++) {
				rt[y+x*r_cols]=0;  // r(y,x)=0;
				for(int z=0;z<cols;z++) 
					// r(y,x)+=(*this)(z,x)*a(y,z);
					rt[y+x*r_cols]+=selft[z+x*self_cols]*at[y+z*a_cols];
			}
		return r;
	}

	MATRIX &operator=(MATRIX r) {
		if(t!=0) delete[] t;
		rows=r.rows;
		cols=r.cols;
		t=new T[rows*cols];
		memcpy(t,r.t,sizeof(T)*rows*cols);
	}

	print() {
		for(int i=0;i<rows;i++) {
			cout << "| ";
			for(int j=0;j<cols;j++) 
				cout << (*this)(j,i) << " ";
			cout << "|\n";
		}
		cout<<"\n";
	}
};


extern "C" { int atoi(const char *); }
void test() {
	MATRIX<float> a(30,50),b(50,30),c;
	SIMPLE_RND rnd;
	for(int x=0;x<30;x++)
		for(int y=0;y<50;y++) {
			a(x,y)=rnd.floating();
			b(y,x)=rnd.floating();
		}

	c=a*b;
	for(int y=0;y<c.rows;y++) {
		for(int x=0;x<c.cols;x++) 
			cout << form("%5d ",(int)(c(x,y)*100.0));
		cout << '\n';
	}
}

void usage() {
	cout << "USAGE: matrix TEST\n"
	        "       matrix [-p] size steps\n";
}

int main(int argc,char *argv[]) {
	SIMPLE_RND rnd;
	int print=0;
	int size=0;
	int steps=0;
	if(argc==2 && strcmp(argv[1],"TEST")==0) {
		test();
		return 0;
	} else if(argc==4 && strcmp(argv[1],"-p")==0) {
		print=1;
		size=atoi(argv[2]);
		steps=atoi(argv[3]);
	} else if(argc==3) {
		size=atoi(argv[1]);
		steps=atoi(argv[2]);
	} else {
		usage();
		return 1;
	}

	MATRIX<float> a(size,size),b(size,size),c;
	for(int i=0;i<size;i++)
		for(int j=0;j<size;j++) {
			a(i,j)=rnd.floating();
			b(i,j)=rnd.floating();
		}
	TIMES t;
	t.start();
	for(int i=0;i<steps;i++) c=a*b;
	TIMES e=t.elapsed();
	if(print) c.print();
	
	cout << e.str();

}