1/*
2 * Benchmark: dgemm / dpotrf / dgetrf through the Fortran-style BLAS/LAPACK
3 * interface. Compiles unchanged against native OpenBLAS and against the
4 * WASM libflame build (f2c built-in BLAS + lapack2flame).
5 *
6 * native: gcc -O2 bench.c -o bench_native -lopenblas
7 * wasm: emcc -O2 bench.c ../install/lib/libflame.a -sALLOW_MEMORY_GROWTH -o bench.js
8 */
14/* The WASM libflame defines these with f2c's int return; native OpenBLAS
15 (gfortran) uses void. Signatures must match exactly under wasm. */
22extern BLASRET dgemm_( char* transa, char* transb, int* m, int* n, int* k,
23 double* alpha, double* a, int* lda, double* b, int* ldb,
24 double* beta, double* c, int* ldc );
25extern int dpotrf_( char* uplo, int* n, double* a, int* lda, int* info );
26extern int dgetrf_( int* m, int* n, double* a, int* lda, int* ipiv, int* info );
28static double now_sec( void )
29{
30 struct timespec ts;
31 clock_gettime( CLOCK_MONOTONIC, &ts );
32 return ts.tv_sec + 1e-9 * ts.tv_nsec;
33}
35/* Deterministic pseudo-random fill so both builds do identical work. */
36static unsigned long long rng_state = 12345;
37static double frand( void )
38{
39 rng_state = rng_state * 6364136223846793005ULL + 1442695040888963407ULL;
40 return ( ( rng_state >> 33 ) & 0xffffff ) / (double) 0x1000000 - 0.5;
41}
43static void fill_random( double* a, int n2 )
44{
45 int i;
46 rng_state = 12345;
47 for ( i = 0; i < n2; i++ ) a[i] = frand();
48}
50/* Diagonally dominant SPD matrix for Cholesky. */
51static void fill_spd( double* a, int n )
52{
53 int i, j;
54 rng_state = 12345;
55 for ( j = 0; j < n; j++ )
56 for ( i = 0; i < n; i++ )
57 a[ j*n + i ] = ( i == j ) ? n : 0.5 * frand();
58 for ( j = 0; j < n; j++ )
59 for ( i = 0; i < j; i++ ) {
60 double v = 0.5 * ( a[ j*n + i ] + a[ i*n + j ] );
61 a[ j*n + i ] = a[ i*n + j ] = v;
62 }
63}
65static void bench_dgemm( int n )
66{
67 double *a = malloc( (size_t)n*n*sizeof(double) );
68 double *b = malloc( (size_t)n*n*sizeof(double) );
69 double *c = malloc( (size_t)n*n*sizeof(double) );
70 double alpha = 1.0, beta = 0.0, t, gflops;
71 char tr = 'N';
73 fill_random( a, n*n );
74 fill_random( b, n*n );
75 memset( c, 0, (size_t)n*n*sizeof(double) );
77 t = now_sec();
78 dgemm_( &tr, &tr, &n, &n, &n, &alpha, a, &n, b, &n, &beta, c, &n );
79 t = now_sec() - t;
81 gflops = 2.0 * n * (double)n * n / t / 1e9;
82 printf( "dgemm n=%5d %10.3f s %8.2f GFLOP/s (check c[0]=%.6f)\n",
83 n, t, gflops, c[0] );
84 free( a ); free( b ); free( c );
85}
87static void bench_dpotrf( int n )
88{
89 double *a = malloc( (size_t)n*n*sizeof(double) );
90 double t, gflops;
91 char lo = 'L';
92 int info = 0;
94 fill_spd( a, n );
95 t = now_sec();
96 dpotrf_( &lo, &n, a, &n, &info );
97 t = now_sec() - t;
99 gflops = ( (double)n * n * n / 3.0 ) / t / 1e9;
100 printf( "dpotrf n=%5d %10.3f s %8.2f GFLOP/s (info=%d, L00=%.6f)\n",
101 n, t, gflops, info, a[0] );
102 free( a );
103}
105static void bench_dgetrf( int n )
106{
107 double *a = malloc( (size_t)n*n*sizeof(double) );
108 int *ipiv = malloc( (size_t)n*sizeof(int) );
109 double t, gflops;
110 int info = 0;
112 fill_random( a, n*n );
113 t = now_sec();
114 dgetrf_( &n, &n, a, &n, ipiv, &info );
115 t = now_sec() - t;
117 gflops = ( 2.0 * n * (double)n * n / 3.0 ) / t / 1e9;
118 printf( "dgetrf n=%5d %10.3f s %8.2f GFLOP/s (info=%d)\n",
119 n, t, gflops, info );
120 free( a ); free( ipiv );
121}
124/* dim_t is int32 in this BLIS build (--int-size=32). */
125extern void bli_thread_set_num_threads( int n_threads );
128int main( int argc, char** argv )
129{
130 int sizes_default[] = { 500, 1000, 2000 };
131 int *sizes = sizes_default, nsizes = 3, i;
133 /* Optional first arg "tN" sets the BLIS thread count (BLIS_MT builds). */
134 if ( argc > 1 && argv[1][0] == 't' ) {
136 int nt = atoi( argv[1] + 1 );
137 if ( nt > 0 ) bli_thread_set_num_threads( nt );
138 printf( "BLIS threads: %d\n", nt );
140 argc--; argv++;
141 }
143 if ( argc > 1 ) {
144 nsizes = argc - 1;
145 sizes = malloc( nsizes * sizeof(int) );
146 for ( i = 0; i < nsizes; i++ ) sizes[i] = atoi( argv[i+1] );
147 }
149 for ( i = 0; i < nsizes; i++ ) bench_dgemm( sizes[i] );
150 for ( i = 0; i < nsizes; i++ ) bench_dpotrf( sizes[i] );
151 for ( i = 0; i < nsizes; i++ ) bench_dgetrf( sizes[i] );
153 return 0;
154}