#define _ISOC11_SOURCE 1
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <malloc.h>

/* aligned_alloc not yet present on g0 */
#define aligned_alloc memalign

extern void matmul(double a[], double b[], double c[], size_t m, size_t n, size_t p);
/* a has m cols and n rows, b has p cols and m rows, c has p cols and n rows */

void init(double a[], size_t n, double start)
{
  size_t i;
  for (i=0; i<n; i++) {
    a[i] = start;
    start += 1.0;
  }
}

void printmatrix(double a[], size_t ncols, size_t nrows)
{
  size_t i,j;
  for (j=0; j<nrows; j++) {
    for (i=0; i<ncols; i++) {
      printf("%7.3g ", a[j*ncols+i]);
    }
    printf("\n");
  }
}

int main(int argc, char *argv[])
{
  if (argc!=2) {
    fprintf(stderr,"Usage: %s n\n",argv[0]);
    exit(1);
  }
  size_t d = atol(argv[1]);
  if (d%4 != 0) {
    fprintf(stderr,"dimension %ld not divisible by 4\n",d);
    exit(1);
  }
  size_t m = d;
  size_t n = d;
  size_t p = d;
  // size_t m = 2;
  // size_t n = 3;
  // size_t p = 4;
  double *a=aligned_alloc(32, m*n*sizeof(double));
  double *b=aligned_alloc(32, m*p*sizeof(double));
  double *c=aligned_alloc(32, n*p*sizeof(double));
  FILE *f = fopen("matmul.out","wb");
  init(a,m*n,1.0);
  // init(b,p*m,8.0);
  init(b,p*m,50000.0);
  // matmul(a,b,c,m,n,p);
  cblas_dgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans,
               n, p, m, 1.0, a, n, b, p, 0.0, c, n);
  // printmatrix(c,p,n);
  fwrite(c, sizeof(double), n*p, f);
  fclose(f);
}