/*
 * Benchmark: dgemm / dpotrf / dgetrf through the Fortran-style BLAS/LAPACK
 * interface. Compiles unchanged against native OpenBLAS and against the
 * WASM libflame build (f2c built-in BLAS + lapack2flame).
 *
 *   native: gcc -O2 bench.c -o bench_native -lopenblas
 *   wasm:   emcc -O2 bench.c ../install/lib/libflame.a -sALLOW_MEMORY_GROWTH -o bench.js
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

/* The WASM libflame defines these with f2c's int return; native OpenBLAS
   (gfortran) uses void. Signatures must match exactly under wasm. */
#ifdef __EMSCRIPTEN__
#define BLASRET int
#else
#define BLASRET void
#endif

extern BLASRET dgemm_( char* transa, char* transb, int* m, int* n, int* k,
                       double* alpha, double* a, int* lda, double* b, int* ldb,
                       double* beta, double* c, int* ldc );
extern int dpotrf_( char* uplo, int* n, double* a, int* lda, int* info );
extern int dgetrf_( int* m, int* n, double* a, int* lda, int* ipiv, int* info );

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

/* Deterministic pseudo-random fill so both builds do identical work. */
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 )
{
  int i;
  rng_state = 12345;
  for ( i = 0; i < n2; i++ ) a[i] = frand();
}

/* Diagonally dominant SPD matrix for Cholesky. */
static void fill_spd( double* a, int n )
{
  int i, j;
  rng_state = 12345;
  for ( j = 0; j < n; j++ )
    for ( i = 0; i < n; i++ )
      a[ j*n + i ] = ( i == j ) ? n : 0.5 * frand();
  for ( j = 0; j < n; j++ )
    for ( i = 0; i < j; i++ ) {
      double v = 0.5 * ( a[ j*n + i ] + a[ i*n + j ] );
      a[ j*n + i ] = a[ i*n + j ] = v;
    }
}

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.3f s  %8.2f GFLOP/s  (check c[0]=%.6f)\n",
          n, t, gflops, c[0] );
  free( a ); free( b ); free( c );
}

static void bench_dpotrf( int n )
{
  double *a = malloc( (size_t)n*n*sizeof(double) );
  double t, gflops;
  char lo = 'L';
  int info = 0;

  fill_spd( a, n );
  t = now_sec();
  dpotrf_( &lo, &n, a, &n, &info );
  t = now_sec() - t;

  gflops = ( (double)n * n * n / 3.0 ) / t / 1e9;
  printf( "dpotrf  n=%5d  %10.3f s  %8.2f GFLOP/s  (info=%d, L00=%.6f)\n",
          n, t, gflops, info, a[0] );
  free( a );
}

static void bench_dgetrf( int n )
{
  double *a = malloc( (size_t)n*n*sizeof(double) );
  int *ipiv  = malloc( (size_t)n*sizeof(int) );
  double t, gflops;
  int info = 0;

  fill_random( a, n*n );
  t = now_sec();
  dgetrf_( &n, &n, a, &n, ipiv, &info );
  t = now_sec() - t;

  gflops = ( 2.0 * n * (double)n * n / 3.0 ) / t / 1e9;
  printf( "dgetrf  n=%5d  %10.3f s  %8.2f GFLOP/s  (info=%d)\n",
          n, t, gflops, info );
  free( a ); free( ipiv );
}

#ifdef BLIS_MT
/* dim_t is int32 in this BLIS build (--int-size=32). */
extern void bli_thread_set_num_threads( int n_threads );
#endif

int main( int argc, char** argv )
{
  int sizes_default[] = { 500, 1000, 2000 };
  int *sizes = sizes_default, nsizes = 3, i;

  /* Optional first arg "tN" sets the BLIS thread count (BLIS_MT builds). */
  if ( argc > 1 && argv[1][0] == 't' ) {
#ifdef BLIS_MT
    int nt = atoi( argv[1] + 1 );
    if ( nt > 0 ) bli_thread_set_num_threads( nt );
    printf( "BLIS threads: %d\n", nt );
#endif
    argc--; argv++;
  }

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

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

  return 0;
}