/*
  y=Ax
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define M 8192 
#define N 8192 
#define L 10

int main()
{
  	int i, j, k;
  	double t, *pa, *px, *py;
 
  	pa = (double *)malloc(M*N*sizeof(double));
  	px = (double *)malloc(N*sizeof(double));
  	py = (double *)malloc(M*sizeof(double));
  	if (!pa||!px||!py) 
    {  printf("Out of memory!\n"); exit(1);  }

	/* set initial value */
	int MN = M*N;
	srand(time(0));
  	for(i=0; i<MN; i++) pa[i] = (double)rand()/RAND_MAX;
  	for(i=0; i<N; i++) px[i] = (double)rand()/RAND_MAX;
  	for(i=0; i<M; i++) py[i] = 0.0;
                
    // i-j loop	    
  	t = (double)clock()/CLOCKS_PER_SEC;
  	for(k=0; k<L; k++)  	
  	for(i=0; i<M; i++) 
  	    for(j=0; j<N; j++)
  	        py[i] += pa[i*N+j] * px[j];  	        
  	t = (double)clock()/CLOCKS_PER_SEC - t;  	
  	printf("row-wise: y[0]=%.2e, y[end]=%.2e, elapsed time: %.2e\n", *py, *(py+M-1), t);

    // j-i loop	    
  	for(i=0; i<M; i++) *(py+i) = 0.0;
  	t = (double)clock()/CLOCKS_PER_SEC;
  	for(k=0; k<L; k++)  	
  	for(j=0; j<N; j++) 
  	    for(i=0; i<M; i++)
  	        py[i] += pa[i*N+j] * px[j];  	        
  	t = (double)clock()/CLOCKS_PER_SEC - t;  	
  	printf("column-wise: y[0]=%.2e, y[end]=%.2e, elapsed time: %.2e\n", *py, *(py+M-1), t);

  	free(pa); free(px); free(py);  
  	return 0;
}