/*
 * Native dgemm benchmark, for a reference point outside the browser.
 * Same flop-count and timing convention as the browser benchmarks.
 *
 *   ./build.sh
 *   OPENBLAS_NUM_THREADS=1 ./bench_native 128 256 512 1024 2048
 *   OPENBLAS_NUM_THREADS=$(nproc) ./bench_native 128 256 512 1024 2048
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

extern void dgemm_( const char* transa, const char* transb, const int* m,
                     const int* n, const int* k, const double* alpha,
                     const double* a, const int* lda, const double* b,
                     const int* ldb, const double* beta, double* c,
                     const int* ldc );

static double now_sec( void )
{
  struct timespec ts;
  clock_gettime( CLOCK_MONOTONIC, &ts );
  return ts.tv_sec + 1e-9 * ts.tv_nsec;
}

/* Same deterministic PRNG as the browser methods, for an identical checksum. */
static unsigned long long rng_state = 12345;
static double frand( void )
{
  rng_state = rng_state * 6364136223846793005ULL + 1442695040888963407ULL;
  return ( ( rng_state >> 33 ) & 0xffffff ) / (double) 0x1000000 - 0.5;
}

static void fill_random( double* a, int n2 )
{
  rng_state = 12345;
  for ( int i = 0; i < n2; i++ ) a[i] = frand();
}

static void bench_dgemm( int n )
{
  double *a = malloc( (size_t)n * n * sizeof(double) );
  double *b = malloc( (size_t)n * n * sizeof(double) );
  double *c = malloc( (size_t)n * n * sizeof(double) );
  double alpha = 1.0, beta = 0.0, t, gflops;
  char tr = 'N';

  fill_random( a, n * n );
  fill_random( b, n * n );
  memset( c, 0, (size_t)n * n * sizeof(double) );

  t = now_sec();
  dgemm_( &tr, &tr, &n, &n, &n, &alpha, a, &n, b, &n, &beta, c, &n );
  t = now_sec() - t;

  gflops = 2.0 * n * (double) n * n / t / 1e9;
  printf( "dgemm  n=%5d  %10.4f s  %8.2f GFLOP/s  (check c[0]=%.6f)\n",
          n, t, gflops, c[0] );

  free( a ); free( b ); free( c );
}

int main( int argc, char** argv )
{
  int sizes_default[] = { 128, 256, 512, 1024, 2048 };
  int *sizes = sizes_default, nsizes = 5;

  if ( argc > 1 ) {
    nsizes = argc - 1;
    sizes = malloc( nsizes * sizeof(int) );
    for ( int i = 0; i < nsizes; i++ ) sizes[i] = atoi( argv[i + 1] );
  }

  for ( int i = 0; i < nsizes; i++ ) bench_dgemm( sizes[i] );

  return 0;
}