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

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

  //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
  for(i=0; i < m; i++)
    for(j=0; j <= n; j++)
      A[i][j] = i + 2*j;
  for(i=0; i < k; i++)
    for(j=0; j <= n; j++)
      B[i][j] = i - 2*j;
  for(i=0; i < m; i++)
    for(j=0; j <= n; j++)
      C[i][j] = 0;

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

  //Print the results
  printf("A is:\n");
  print_mat(A,m,k);
  printf("B is:\n");
  print_mat(B,k,n);
  printf("C is:\n");
  print_mat(C,m,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");
  }
}