/*
  This matrix multiply driver was originally written by Jason Riedy.

  The output format: "Size: %u\tmflop/s: %g\n"
*/


#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <float.h>
#include <math.h>

#include <sys/types.h>
#include <sys/resource.h>

#include <unistd.h>
#include <time.h>

#define MIN_RUNS     4
#define MIN_CPU_SECS 1.0
#define MAX_SIZE 512

const unsigned test_sizes[] = {
     20,
     24,
     31,
     32,
     48,
     64,
     73,
     96,
     97,
     127,
     128,
     129,
     163,
     191,
     192,
     229,
     255,
     256,
     257,
     319,
     320,
     321,
     417,
     479,
     480,
     511,
     512,
};


#define N_SIZES ((unsigned) sizeof (test_sizes) / sizeof (unsigned))

double A[MAX_SIZE][MAX_SIZE];
double B[MAX_SIZE][MAX_SIZE];
double C[MAX_SIZE][MAX_SIZE];


void matrix_init (double A[MAX_SIZE][MAX_SIZE])
{
     unsigned row, col;

     for (row = 0; row < MAX_SIZE; row++) {
         for (col = 0; col < MAX_SIZE; col++) {
            A[row][col] = drand48 ();
	 }
     }
}

void matrix_clear (double C[][]) 
{
     memset (C, 0, MAX_SIZE * MAX_SIZE * sizeof (double));
}

void mult_mat(const unsigned M,
              double A[MAX_SIZE][MAX_SIZE], 
	      double B[MAX_SIZE][MAX_SIZE], 
	      double C[MAX_SIZE][MAX_SIZE])
{
    unsigned i, j, k;

    matrix_clear (C);

    for (i = 0; i < M; ++i) {
        C[i][i] = A[i][i] + B[i][i];
    }
}

double timer (const unsigned M,
              double A[MAX_SIZE][MAX_SIZE], 
	      double B[MAX_SIZE][MAX_SIZE], 
	      double C[MAX_SIZE][MAX_SIZE])
{

    clock_t cpu_time;
    clock_t last_clock;
    double mflops, mflop_s;
    double secs = -1;

    unsigned num_iterations = MIN_RUNS;
    unsigned i;

    cpu_time = 0;

    matrix_clear (C);
    last_clock = clock();
    for (i = 0; i < num_iterations; ++i) {
         mult_mat(M, A, B, C);
    }
    cpu_time += clock() - last_clock;

//    mflops  = 2.0 * num_iterations * M * M * M / 1.0e6;
    mflops  = 2.0 * num_iterations * M / 1.0e6;
    secs    = cpu_time / ((double) CLOCKS_PER_SEC);
    mflop_s = mflops/secs;

printf("Secs: %lf\t", secs);

    return mflop_s;
}

int main ()
{
     unsigned sz_i;
     double mflop_s;

     matrix_init (A);
     matrix_init (B);

     for (sz_i = 0; sz_i < N_SIZES; ++sz_i) {

          const unsigned M = test_sizes[sz_i];

          mflop_s = timer(M, A, B, C);    

          printf ("Size: %u\tmflop/s: %g\n", M, mflop_s);
     }

     return 0;
}