#include <mach/mach.h>
#include <mach/mach_time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "mkl.h"

const int matrix_size[] = {
     10,    13,    16,    21,    26,    34,    43,    55,
     70,    89,   113,   144,   183,   234,   298,   379,
    483,   616,   785,  1000
};

const int matrix_size_pow2[] = {
    4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048
};

const int NREPEATS = 30;
const int MAX_N = 2048;

double nanodiff(const uint64_t _t0, const uint64_t _t1, 
                  const mach_timebase_info_data_t *tb_info)
{   
    long double t0, t1, numer, denom, nanosec;
    numer = (long double)(tb_info->numer);
    denom = (long double)(tb_info->denom);    
    t0 = (long double)(_t0);
    t1 = (long double)(_t1);
    nanosec = (t1 - t0) * numer / denom;
    return (double)nanosec;
}

void fill_with_random(const int n, double *x)
{
     static unsigned int m_w = 123456;
     static unsigned int m_z = 5635273;
     int i;
     for (i=0; i<n; i++) {
         m_z = 36969 * (m_z & 65535) + (m_z >> 16);
         m_w = 18000 * (m_w & 65535) + (m_w >> 16);
         *x++ = ((m_z << 16) + m_w) * 2.3283064365386963e-10;
     }
}


void statistics(const int n, const double *x, 
                   double *m, double *s, 
                   double *min, double *max)
{
    double sum_x=0.0, cx=0.0, sum_cxcx=0.0, _m;
    double minval, maxval, v;
    int i;
    for (i=0; i<n; i++) sum_x += x[i];
    _m = sum_x / (double)n;
    for (i=0; i<n; i++) {
        cx = x[i] - _m;
        sum_cxcx += cx*cx;
    }    
    *m = _m;
    *s = sqrt(sum_cxcx / (double)(n-1));
    minval = *x;
    maxval = *x;
    for (i=1; i<n; i++) {
        v = *x++;
        maxval = (maxval < v ? v : maxval);
        minval = (minval > v ? v : minval);
    }
    *max = maxval;
    *min = minval;
}


int main(int argc, char **argv)
{

    double nanosec[NREPEATS];
    uint64_t t0, t1;
    mach_timebase_info_data_t tb_info;
    
    double *A = (double*)mkl_malloc(MAX_N*MAX_N*sizeof(double),64); 
    double *B = (double*)mkl_malloc(MAX_N*MAX_N*sizeof(double),64); 
    double *C = (double*)mkl_malloc(MAX_N*MAX_N*sizeof(double),64);
    double mean, std, min, max;
    
    int i, j, k, m, n;
    
    mach_timebase_info(&tb_info);
    
    fill_with_random(MAX_N*MAX_N, A);
    fill_with_random(MAX_N*MAX_N, B);
    fill_with_random(MAX_N*MAX_N, C);
        
    for (i=0; i<20; i++) {
        n = matrix_size[i];
        m = n; k = n;
        for (j=0; j<NREPEATS; j++) {
    
            t0 = mach_absolute_time();
    
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                m, n, k, 1.0, A, k, B, n, 1.0, C, n);

            t1 = mach_absolute_time();
    
            nanosec[j] = nanodiff(t0, t1, &tb_info);
        }
        
        statistics(NREPEATS, nanosec, &mean, &std, &min, &max);
        
        printf("[%4d,    %.4e,    %.4e,    %.4e,    %.4e],\n", 
                n, mean, std, min, max);
        
    }
    
    printf("\n\n");
    
    for (i=0; i<10; i++) {
        n = matrix_size_pow2[i];
        m = n; k = n;
        for (j=0; j<NREPEATS; j++) {
    
            t0 = mach_absolute_time();
    
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                m, n, k, 1.0, A, k, B, n, 1.0, C, n);

            t1 = mach_absolute_time();
    
            nanosec[j] = nanodiff(t0, t1, &tb_info);
        }
        
        statistics(NREPEATS, nanosec, &mean, &std, &min, &max);
        
        printf("[%4d,    %.4e,    %.4e,    %.4e,    %.4e],\n", 
                n, mean, std, min, max);        
    }
        
    mkl_free(A); mkl_free(B); mkl_free(C);
}