#include <stdlib.h>         // for malloc(), free(), rand(), srand()
#include <stdio.h>          // for printf()
#include <stddef.h>         // for size_t, ptrdiff_t
#include <math.h>           // for nan(), fabs()
#include <float.h>          // for DBL_EPSILON
#include <stdbool.h>        // for typedef bool
#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
randGeMatrix(size_t m, size_t n, bool withNan,
             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] = withNan
                                     ? nan("")
                                     : 2.*(rand()-RAND_MAX/2)/RAND_MAX;
        }
    }
}

void
printGeMatrix(size_t m, size_t n,
              const double *A,
              ptrdiff_t incRowA, ptrdiff_t incColA)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            printf("%9.2lf ", A[i*incRowA + j*incColA]);
        }
        printf("\n");
    }
    printf("\n");
}

double
dgenrm_inf(size_t m, size_t n,
           const double *A,
           ptrdiff_t incRowA, ptrdiff_t incColA)
{
    double result = 0;

    for (size_t i=0; i<m; ++i) {
        double sum = 0;

        for (size_t j=0; j<n; ++j) {
            sum += fabs(A[i*incRowA+j*incColA]);
        }

        if (sum>result) {
            result = sum;
        }
    }

    return result;
}

//-- BLAS Level 1 functions ----------------------------------------------------

void
dcopy(size_t n,
      const double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    for (size_t i=0; i<n; ++i) {
        y[i*incY] = x[i*incX];
    }
}

void
daxpy(size_t n, double alpha,
      const double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    for (size_t i=0; i<n; ++i) {
        y[i*incY] += alpha*x[i*incX];
    }
}

double
ddot(size_t n,
     const double *x, ptrdiff_t incX,
     const double *y, ptrdiff_t incY)
{
    double result = 0;

    for (size_t i=0; i<n; ++i) {
        result += x[i*incX]*y[i*incY];
    }

    return result;
}

void
dscal(size_t n,
      double alpha,
      double *x, ptrdiff_t incX)
{
    if (alpha==1) {
        return;
    }
    if (alpha!=0) {
        for (size_t i=0; i<n; ++i) {
            x[i*incX] *= alpha;
        }
    } else {
        for (size_t i=0; i<n; ++i) {
            x[i*incX] = 0;
        }
    }
}

//-- BLAS Level 2 functions ----------------------------------------------------

void
dgemv_dot(size_t m, size_t n,
          double alpha,
          const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
          const double *x, ptrdiff_t incX,
          double beta,
          double *y, ptrdiff_t incY)
{
    /* Your code here */
}

void
dgemv_axpy(size_t m, size_t n,
           double alpha,
           const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
           const double *x, ptrdiff_t incX,
           double beta,
           double *y, ptrdiff_t incY)
{
    /* your code here */
}


#ifndef DOTF
#define DOTF 2
#endif

void
dgemv_dotf(size_t m, size_t n,
           double alpha,
           const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
           const double *x, ptrdiff_t incX,
           double beta,
           double *y, ptrdiff_t incY)
{
    /* your code here */
}


#ifndef AXPYF
#define AXPYF 2
#endif

void
dgemv_axpyf(size_t m, size_t n,
            double alpha,
            const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
            const double *x, ptrdiff_t incX,
            double beta,
            double *y, ptrdiff_t incY)
{
    /* Your code here */
}


//-- BLAS Level 2: dgemv reference implementation and error bound --------------

void
dgemv_ref(size_t m, size_t n,
          double alpha,
          const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
          const double *x, ptrdiff_t incX,
          double beta,
          double *y, ptrdiff_t incY)
{
    if (beta!=1) {
        if (beta!=0) {
            for (size_t i=0; i<m; ++i) {
                y[i*incY] *= beta;
            }
        } else {
            for (size_t i=0; i<m; ++i) {
                y[i*incY] = 0;
            }
        }
    }
    if (alpha!=0) {
        for (size_t j=0; j<n; ++j) {
            for (size_t i=0; i<m; ++i) {
                y[i*incY] += alpha*A[i*incRowA+j*incColA]*x[j*incX];
            }
        }
    }
}

// - Computes error bound for the test result ySol of the gemv operation
//   beta*y0 + alpha*A*x.
// - yRef is trusted result.
// - ySol gets overwritten.
double
dgemv_err(size_t m, size_t n,
          double alpha,
          const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
          const double *x, ptrdiff_t incX,
          const double *y0, ptrdiff_t incY0,
          double beta,
          const double *yRef, ptrdiff_t incYRef,
          double *ySol, ptrdiff_t incYSol)
{
    double nrmY0 = dgenrm_inf(m, 1, y0, incY0, 1);
    double nrmX  = dgenrm_inf(n, 1, x, incX, 1);
    double nrmA  = dgenrm_inf(m, n, A, incRowA, incColA);
    size_t maxMN = m<n ? n : m;

    // nrmDiff = ||y2 - y1||_inf
    daxpy(m, -1, yRef, incYRef, ySol, incYSol);
    double nrmDiff = dgenrm_inf(m, 1, ySol, incYSol, 1);

    return nrmDiff /
        (DBL_EPSILON*(maxMN*fabs(alpha)*nrmA*nrmX + m*fabs(beta)*nrmY0));
}



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

#ifndef COLMAJOR
#define COLMAJOR 0
#endif

#ifndef SEED_RAND
#define SEED_RAND 0
#endif

#ifndef MAX_M
#define MAX_M 4000
#endif

#ifndef MAX_N
#define MAX_N 4000
#endif

#ifndef ALPHA
#define ALPHA 1
#endif

#ifndef BETA
#define BETA 1
#endif

int
main()
{
    srand(SEED_RAND);
    printf("#COLMAJOR = %d\n", COLMAJOR);
    printf("#ALPHA    = %lf\n", (double)ALPHA);
    printf("#BETA     = %lf\n", (double)BETA);

    double *A    = malloc(MAX_M*MAX_N*sizeof(double));
    double *x    = malloc(MAX_N*sizeof(double));
    double *y0   = malloc(MAX_M*sizeof(double));
    double *yRef = malloc(MAX_M*sizeof(double));
    double *ySol = malloc(MAX_M*sizeof(double));
    if (!A || !x || !y0 || !yRef || !ySol) {
        abort();
    }

    // print header
    printf("#%4s %4s ", "m", "n");
    printf("%10s %10s ", "time ref", "mflops ref");
    printf("%10s %10s %7s ", "time 1", "mflops 1", "err");
    printf("%10s %10s %7s ", "time 2", "mflops 2", "err");
    printf("\n");

    for (size_t m=10, n=10; m<=MAX_M && n<=MAX_N; m+=10, n+=10) {
        ptrdiff_t incRowA = COLMAJOR ? 1 : n;
        ptrdiff_t incColA = COLMAJOR ? m : 1;

        ptrdiff_t incX    = 1;
        ptrdiff_t incY0   = 1;
        ptrdiff_t incYRef = 1;
        ptrdiff_t incYSol = 1;

        double    alpha = ALPHA;
        double    beta  = BETA;

        double    mflop = m*(2*n+1)/1000./1000.;

        randGeMatrix(m, n, ALPHA==0, A, incRowA, incColA);
        randGeMatrix(n, 1, ALPHA==0, x, incX, 0);
        randGeMatrix(m, 1, BETA==0, y0, incY0, 0);

        printf("A =\n");
        printGeMatrix(m, n, A, incRowA, incColA);

        printf("x =\n");
        printGeMatrix(1, n, x, 0, incX);

        printf("y0 =\n");
        printGeMatrix(1, m, y0, 0, incY0);

        printf(" %4zu %4zu ", m,  n);

        {
            double t    = 0;
            size_t runs = 0;

            while (t<0.1 || runs<3) {
                dcopy(m, y0, incY0, yRef, incYRef);

                double t0 = walltime();
                dgemv_ref(m, n,
                          alpha,
                          A, incRowA, incColA,
                          x, incX,
                          beta,
                          yRef, incYRef);
                t += walltime() - t0;
                ++runs;
            }
            t /= runs;

            printf("\nyRef =\n");
            printGeMatrix(1, m, yRef, 0, incYRef);
            printf("%10.2lf %10.2lf ", t, mflop/t);
        }

        {
            double t    = 0;
            size_t runs = 0;

            while (t<0.1 || runs<3) {
                dcopy(m, y0, incY0, ySol, incYSol);

                double t0 = walltime();
                dgemv_dot(m, n,
                          alpha,
                          A, incRowA, incColA,
                          x, incX,
                          beta,
                          ySol, incYSol);
                t += walltime() - t0;
                ++runs;
            }
            t /= runs;

            printf("\nySol =\n");
            printGeMatrix(1, m, ySol, 0, incYSol);

            double err = dgemv_err(m, n, alpha,
                                   A, incRowA, incColA,
                                   x, incX,
                                   y0, incY0,
                                   beta,
                                   yRef, incYRef,
                                   ySol, incYSol);

            printf("\nyDiff =\n");
            printGeMatrix(1, m, ySol, 0, incYSol);

            printf("%10.2lf %10.2lf %7.1e ", t, mflop/t, err);
        }

        printf("\n");
        break;
    }

    free(A);
    free(x);
    free(y0);
    free(yRef);
    free(ySol);
}