#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <cstring>
#include <cmath>
#include <chrono>
#include <ctime>

using namespace std;

//Declare BLAS functions without looking for them 
extern "C" {
  extern void dgemm_(char*, char*, int*, int*, int*, double*, double*, int*, double*, int*, double*, double*, int*);
}


double **matrix_init(int row, int col);
void matrix_delete(double **A);
void print_mat(double **A, int rows, int cols);

int main(int argc, char *argv[])
{
  int m, n, k; //matrix sizes

  if(argc < 2) {
    printf("Must supply matrix size (m or m n k) as argument.\n");
    return 1;
  }
  m = atoi(argv[1]);
  if(argc != 4) {
    n = m;
    k = m;
  } else {
    n = atoi(argv[2]);
    k = atoi(argv[3]);
  }

  int i, j; //iterators

  //timer
  chrono::time_point<chrono::system_clock> start, end;
  
  //inputs to dgemm
  char transa = 'n', transb = 'n';
  double alpha = 1.0, beta = 0.0;

  // Allocate and fill matrices
  double **A = matrix_init(m, k);
  double **B = matrix_init(k, n);
  double **C = matrix_init(m, n);

  // Initialize random seed
  srand(time(NULL));

  // Initialize
  for(i=0; i < m; i++)
    for(j=0; j <= n; j++)
      A[i][j] = ((double) rand())/((double) RAND_MAX);
  for(i=0; i < k; i++)
    for(j=0; j <= n; j++)
      B[i][j] = ((double) rand())/((double) RAND_MAX);
  for(i=0; i < m; i++)
    for(j=0; j <= n; j++)
      C[i][j] = 0;

  //Get system time
  start = chrono::system_clock::now();
  time_t start_time = chrono::system_clock::to_time_t(start);
  cout << "started computation at " << ctime(&start_time);

  //Multiply matrices. Output is returned in c
  dgemm_(&transb, &transa, &n, &m, &k, &alpha, B[0], &n, A[0], &k, &beta, C[0], &n);

  //Get new system time and print difference
  end = chrono::system_clock::now();
  chrono::duration<double> elapsed_seconds = end-start;
  time_t end_time = chrono::system_clock::to_time_t(end);

  cout << "finished computation at " << ctime(&end_time)
            << "elapsed time: " << elapsed_seconds.count() << "s\n";

  cout << "Gflops/s for MM was " << (((2.0*((double) m)*((double) n)*((double) k))/pow(10.0,9))/elapsed_seconds.count()) << "\n";

  //Free the memory
  matrix_delete(A);
  matrix_delete(B);
  matrix_delete(C);

  return 0;
}


// Allocate a (row-wise) matrix of dimension (row x col)
// Rows are contiguous in memory
double **matrix_init(int row, int col)
{
  double **A = new double* [row];
  double *B = new double[row*col];
  memset(static_cast<void*>(B), 0, row*col*sizeof(double));
  for(int i=0; i < row; i++)
    A[i] = &(B[i*col]);
  return A;
}

// Delete a matrix A of dimension dim initialized with matrix_init()
void matrix_delete(double **A)
{
  delete[] A[0];
  delete[] A;
}

// Print a matrix to the screen
void print_mat(double **A, int rows, int cols)
{
  for(int i=0; i<rows; i++)
  {
    for(int j=0; j<cols; j++)
      printf("%8.3f ",A[i][j]);

    printf("\n");
  }
}