public class MatMul{

    public static void main(String args[]) {
	PAPICounter t1 = new PAPICounter(PAPICounter.PAPI_FLOPS);
	PAPICounter t2 = new PAPICounter(PAPICounter.PAPI_FP_INS);
	t1.start();
	t2.start();
	matmulDriver();
	t1.stop();
	t2.stop();
	t2.start();
	matmulDriver();
	t2.stop();
	System.out.println("FLOPS: " + t1.getCounterValue() + " on processor " + 
			   Ti.thisProc());
	System.out.println("Total number of floating point operations: " + 
			   t2.getCounterValue() + " on processor " + Ti.thisProc());
    }
    
    public static void matMul(double [2d] a, double [2d] b, double [2d] c) {

	RectDomain<2> aRow; 
	RectDomain<2> bCol; 
	Point<2> ij;
	Point<2> ik;
	Point<2> kj;
	int n = c.domain().max()[1];  // assumes square for now
	// why are points 1-based?
       
	double sum;
	
	foreach (ij within c.domain()) {
	    int i = ij[1];
	    int j = ij[2];
	    aRow = [i:i,0:n];
	    
	    foreach (ik within aRow) {
		kj = [ik[2],ij[2]];
	        c[ij] += a[ik] * b[kj];
	    }
	}
    }

    public static void matmulDriver () {

	for (int n = 32; n <= 32; n*=2) {
	    
	    RectDomain<2> d2 = [0:n-1,0:n-1];
	    double [2d] a = new double [d2];
	    double [2d] b = new double [d2];
	    double [2d] c = new double [d2];
	    foreach (ij within d2) {
		a[ij] = 2.0;
		b[ij] = 2.0;
		c[ij] = 0.0;
	    }
	    
	    long ms = System.currentTimeMillis();
	    matMul(a,b,c);
	    ms = System.currentTimeMillis() - ms;
	    //System.out.println("n = " + n + ", Mflops = " + 
		//	       ((double) 2*n*n*n)/((double) ms*1000));
	    foreach (ij within d2) {
		if (c[ij] != 4*n) System.out.println("Incorrect computation");
	    }
	}
    }
}