#include <math.h>           // for nan()
#include <stddef.h>         // for size_t, ptrdiff_t
#include <stdio.h>          // for printf
#include <stdlib.h>         // for malloc(), free()
#include <sys/times.h>      // needed for walltime()
#include <unistd.h>         // needed for walltime()

//-- Function for benchmarking and testing -------------------------------------

double
walltime()
{
   struct tms    ts;
   static double ClockTick=0.0;

   if (ClockTick==0.0) {
        ClockTick = 1.0 / ((double) sysconf(_SC_CLK_TCK));
   }
   return ((double) times(&ts)) * ClockTick;
}

//------------------------------------------------------------------------------

void
initGeMatrix(size_t m, size_t n,
             double *A,
             ptrdiff_t incRowA, ptrdiff_t incColA)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            A[i*incRowA + j*incColA] = i*n + j +1;
        }
    }
}

void
dgeaxpy(size_t m, size_t n,
        double alpha,
        const double *A,
        ptrdiff_t incRowA, ptrdiff_t incColA,
        double *B,
        ptrdiff_t incRowB, ptrdiff_t incColB)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            B[i*incRowB + j*incColB] += alpha*A[i*incRowA + j*incColA];
        }
    }
}

#define MAX_M 2*16*16*64
#define MAX_N 2*16*16*64

int
main()
{
    double *A = malloc(MAX_M*MAX_N*sizeof(double));
    double *B = malloc(MAX_M*MAX_N*sizeof(double));

    printf("#%5s %5s %10s %10s\n", "m", "n", "t (colm)", "t (rowm)");

    for (size_t m=128, n=128; m<=MAX_M && n<=MAX_N; m+=16, n+=16) {

        printf(" %5ld %5ld ", (long)m, (long)n);


        // bench initialization for col major matrix
        {

            ptrdiff_t incRowA = MAX_N;
            ptrdiff_t incColA = 1;

            ptrdiff_t incRowB = MAX_N;
            ptrdiff_t incColB = 1;

            double    alpha = 3;

            initGeMatrix(m, n, A, incRowA, incColA);
            double t  = walltime();
            dgeaxpy(m, n, alpha, A, incRowA, incColA, B, incRowB, incColB);
            t = walltime() - t;

            printf("%10.2lf ", t);
        }

        // bench initialization for col major matrix
        {

            ptrdiff_t incRowA = 1;
            ptrdiff_t incColA = MAX_M;

            ptrdiff_t incRowB = 1;
            ptrdiff_t incColB = MAX_M;

            double    alpha = 3;

            initGeMatrix(m, n, A, incRowA, incColA);
            double t  = walltime();
            dgeaxpy(m, n, alpha, A, incRowA, incColA, B, incRowB, incColB);
            t = walltime() - t;

            printf("%10.2lf ", t);
        }

        printf("\n");

    }
    free(B);
    free(A);
}