/*
 * Browser benchmark entry point for libflame+BLIS WASM.
 * Exposes run_bench(routine, n, nthreads) for calling from JavaScript;
 * returns GFLOP/s (or -1 on error) and printf's a detail line (captured
 * by the page via Module.print).
 *
 * Built two ways by build-web.sh:
 *   single-threaded: BLIS st + libflame.a
 *   threaded:        -pthread -DBLIS_MT, BLIS pthreads + libflame-mt.a
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <emscripten.h>

extern int 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 );

#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

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

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();
}

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 double 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';

  if ( !a || !b || !c ) { free(a); free(b); free(c); return -1.0; }
  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  %8.3f s  %7.2f GFLOP/s  (check c[0]=%.6f)\n",
          n, t, gflops, c[0] );
  free( a ); free( b ); free( c );
  return gflops;
}

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

  if ( !a ) return -1.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  %8.3f s  %7.2f GFLOP/s  (info=%d, L00=%.6f)\n",
          n, t, gflops, info, a[0] );
  free( a );
  return info == 0 ? gflops : -1.0;
}

static double 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;

  if ( !a || !ipiv ) { free(a); free(ipiv); return -1.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  %8.3f s  %7.2f GFLOP/s  (info=%d)\n",
          n, t, gflops, info );
  free( a ); free( ipiv );
  return info == 0 ? gflops : -1.0;
}

EMSCRIPTEN_KEEPALIVE
double run_bench( int routine, int n, int nthreads )
{
  if ( n < 2 || n > 8000 ) return -1.0;

#ifdef BLIS_MT
  if ( nthreads > 0 ) bli_thread_set_num_threads( nthreads );
#else
  (void) nthreads;
#endif

  switch ( routine ) {
    case 0: return bench_dgemm( n );
    case 1: return bench_dpotrf( n );
    case 2: return bench_dgetrf( n );
    default: return -1.0;
  }
}