#include <u.h>
#include <libn.h>
#include <libmath.h>
#include <libmath/blas.h>
#include <time.h>

#include <vendor/blas/cblas.h>

#define NCOL 2000
#define NROW 2000
#define NIT  2000
#define INC  2
error
main()
{
    int i, j, nit;
    double *x, *y, *m, res[2];

    clock_t t;
    double tprof[2] = { 0 };

    rng·init(0);

    x = malloc(sizeof(*x)*NCOL);
    y = malloc(sizeof(*x)*NROW);
    m = malloc(sizeof(*x)*NROW*NCOL);

#define DO_0  t = clock();                                  \
              blas·gerd(NROW/INC, NCOL/INC, 1.2, x, INC, y, INC, m, NCOL); \
              res[0] += m[0]; \
              t = clock() - t;                              \
              tprof[0] += 1000.*t/CLOCKS_PER_SEC;           \

#define DO_1 t = clock();                                   \
             cblas_dger(CblasRowMajor, NROW/INC, NCOL/INC, 1.2,  x, INC, y, INC, m, NCOL); \
             res[1] += m[0]; \
             t = clock() - t;                               \
             tprof[1] += 1000.*t/CLOCKS_PER_SEC;      

    for (nit = 0; nit < NIT; nit++) {
        for (i = 0; i < NROW; i++) {
            x[i] = rng·random();
            y[i] = rng·random();
            for (j = 0; j < NCOL; j++) {
                m[j + NCOL*i] = rng·random();
            }
        }

        switch (nit % 2) {
            case 0: DO_0; DO_1; break;
            case 1: DO_1; DO_0; break;
        }
    }
    printf("mean time/iteration (mine): %fms\n", tprof[0]/NIT);
    printf("--> result (mine): %f\n", res[0]);
    printf("mean time/iteration (openblas): %fms\n", tprof[1]/NIT);
    printf("--> result (openblas): %f\n", res[1]);

    return 0;
}