/* matrix multiply permutations */
#include <stdio.h>
#include <sys/resource.h>
#include <sys/time.h>

#define MAXN 1024  
#define MINN 64  
#define INCN 64

/* global float arrayss */
float ga[MAXN][MAXN], gb[MAXN][MAXN], gc[MAXN][MAXN]; 
void ijk(float [][MAXN] , float [][MAXN] , float [][MAXN] , int );
void jik(float [][MAXN] , float [][MAXN] , float [][MAXN] , int );
void jki(float [][MAXN] , float [][MAXN] , float [][MAXN] , int );
void kji(float [][MAXN] , float [][MAXN] , float [][MAXN] , int );
void kij(float [][MAXN] , float [][MAXN] , float [][MAXN] , int );
void ikj(float [][MAXN] , float [][MAXN] , float [][MAXN] , int );

void uswtime(double *usertime, double *systime, double *walltime)
{
double mega = 1.0e-6;
struct rusage buffer;
struct timeval tp;
struct timezone tzp;

getrusage(RUSAGE_SELF, &buffer);
gettimeofday(&tp, &tzp);
*usertime = (double) buffer.ru_utime.tv_sec +
              1.0e-6 * buffer.ru_utime.tv_usec;
*systime  = (double) buffer.ru_stime.tv_sec +
              1.0e-6 * buffer.ru_stime.tv_usec;
*walltime = (double) tp.tv_sec + 1.0e-6 * tp.tv_usec;
}


/* check the result float [][] for correctness */
void checkresult(float c[][MAXN], int n)
{
    int i, j;

    for (i = 0; i < n; i++)
	for (j = 0; j < n; j++)
	    if (c[i][j] != (double)n) {
		printf("Error: bad number (%f) in result matrix (%d,%d)\n", 
		       c[i][j], i, j);
		exit(0);
	    }
}

/* Run f and return clocks per inner loop iteration */
double run(int ftype, int n)
{
double    utime0, stime0, wtime0,
          utime1, stime1, wtime1;

    
     uswtime(&utime0, &stime0, &wtime0);
switch(ftype) {
    case 1:ijk(ga, gb, gc, n);break;
    case 2:jik(ga, gb, gc, n);break;
    case 3:jki(ga, gb, gc, n);break;
    case 4:kji(ga, gb, gc, n);break;
    case 5:kij(ga, gb, gc, n);break;
    case 6:ikj(ga, gb, gc, n);break;
    }	  
     uswtime(&utime1, &stime1, &wtime1);
     return(utime1 - utime0 + stime1 - stime0);
}

/* reset result float [][] to zero */
void reset(float c[][MAXN], int n)
{
    int i,j;
    
    for (i = 0; i < n; i++) {
	for (j = 0; j < n; j++) {
	    c[i][j] = 0.0;
	}
    }
}

/* initialize input float [][]s to 1 */
void init(float a[][MAXN], float b[][MAXN], int n) 
{
    int i,j;
    
    for (i = 0; i < n; i++) {
	for (j = 0; j < n; j++) {
	    a[i][j] = 1.0;
	    b[i][j] = 1.0;
	}
    }
}


/* print an float [][] (debug) */
void printarray(float a[][MAXN], int n)
{
    int i, j;

    for (i = 0; i < n; i++) {
	for (j = 0; j < n; j++) {
	    printf("%5.1f ", a[i][j]);
	}
	printf("\n");
    }
}

/***********************************************
 * Six different versions of matrix multiply 
 ***********************************************/
void ijk(float A[][MAXN], float B[][MAXN], float C[][MAXN], int n)
{
    int i, j, k;
    double sum;

    /* $begin mm-ijk */
for (i = 0; i < n; i++) 
    for (j = 0; j < n; j++) {
	sum = 0.0;
	for (k = 0; k < n; k++)
	    sum += A[i][k]*B[k][j];
	C[i][j] += sum;
    }
/* $end mm-ijk */

}

void jik(float A[][MAXN], float B[][MAXN], float C[][MAXN], int n)
{
    int i, j, k;
    double sum;

/* $begin mm-jik */
for (j = 0; j < n; j++) 
    for (i = 0; i < n; i++) {
	sum = 0.0;
	for (k = 0; k < n; k++)
	    sum += A[i][k]*B[k][j];
	C[i][j] += sum;
    }
/* $end mm-jik */
}

void ikj(float A[][MAXN], float B[][MAXN], float C[][MAXN], int n)
{
    int i, j, k;
    double r;
    
    /* $begin mm-ikj */
for (i = 0; i < n; i++)
    for (k = 0; k < n; k++) {
	r = A[i][k];
	for (j = 0; j < n; j++)
	    C[i][j] += r*B[k][j];
    }
/* $end mm-ikj */
}

void kij(float A[][MAXN], float B[][MAXN], float C[][MAXN], int n)
{
    int i, j, k;
    double r;

    /* $begin mm-kij */
for (k = 0; k < n; k++)
    for (i = 0; i < n; i++) {
	r = A[i][k];
	for (j = 0; j < n; j++)
	    C[i][j] += r*B[k][j];
    }
/* $end mm-kij */
}

void kji(float A[][MAXN], float B[][MAXN], float C[][MAXN], int n)
{
    int i, j, k;
    double r;

/* $begin mm-kji */
for (k = 0; k < n; k++)
    for (j = 0; j < n; j++) {
	r = B[k][j];
	for (i = 0; i < n; i++)
	    C[i][j] += A[i][k]*r;
    }
/* $end mm-kji */
}

void jki(float A[][MAXN], float B[][MAXN], float C[][MAXN], int n)
{
    int i, j, k;
    double r;

/* $begin mm-jki */
for (j = 0; j < n; j++)
    for (k = 0; k < n; k++) {
	r = B[k][j];
	for (i = 0; i < n; i++)
	    C[i][j] += A[i][k]*r;
    }
/* $end mm-jki */
}

/* 
 * Run the six versions of matrix multiply and display performance
 * as clock cycles per inner loop iteration.
 */ 

int main()
{
    int n;

    init(ga, gb, MAXN);

    printf("Tiempo de Computacion\n");
    printf("%3s%6s%6s%6s%6s%6s%6s\n", "n", 
	   "ijk", "jik", "jki", "kji", "kij", "ikj");

    for (n = MINN; n <= MAXN; n += INCN) {  
	printf("%3d ", n);

	printf("%5.2f ", run(1, n));
	printf("%5.2f ", run(2, n));
	printf("%5.2f ", run(3, n));
	printf("%5.2f ", run(4, n));
	printf("%5.2f ", run(5, n));
	printf("%5.2f ", run(6, n));
	printf("\n");
    }
    exit(0);
}