/* matrix_multiply.c
 * 
 *   $ gcc -O2 matrix_multiply.c -o matrix_multiply
 *   $ time ./matrix_multiply
 *   The average C value is 12645.7
 *   real	2m29.296s
 *   user	2m29.263s
 *   sys	0m0.004s
 *
 * Computer Systems course, Nov 2020, cs.bennington.college
 */

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

#define N 1000
#define REPEATS 100
#define randint(range) ((rand() % (2*(range))) - (range))

float A[N][N];
float B[N][N];
float C[N][N];

int main(){
  int i, j, k, r;
  float sum, total;

  total = 0;
  for (r=0; r<REPEATS; r++){
  
    // First put random values into A and B. 
    for (i=0; i<N; i++){ for (j=0; j<N; j++) A[i][j] = 1.0*randint(16); }
    for (i=0; i<N; i++){ for (j=0; j<N; j++) B[i][j] = 1.0*randint(32); }
    
    // Now multiply them together using one one of the methods
    // given in Figure 6.44 of our Computer Systems text.
    //
    // Replace this part with the other approaches in Figure 6.44,
    // compare runtimes, and discuss in terms of locality, memory,
    // processor speed, and all that. 
    //
    // ----------- version 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;                                     // --
      }                                                     // --
    }                                                       // --
    // ----------------------------------------------------------
    
    // Just for the heck of it, accumulate the sum of C's values.
    sum = 0.0;
    for (i=0; i<N; i++){ for (j=0; j<N; j++) sum += C[i][j]; }
    total += sum / (N*N);
  }
  printf("The average C value is %g \n", total/REPEATS);
  return 0;
}