#include <stdio.h>
#include <bench/aux.h>
#include <bench/refblas.h>
#include <bench/errbound.h>
#include <ulmblas/level1.h>
#include <ulmblas/level2.h>
#include <ulmblas/level3.h>

//-- TRLSM with simple cache optimization --------------------------------------

#ifndef DTRLSM_MC
#define DTRLSM_MC 64
#endif

#ifndef DTRLSM_NC
#define DTRLSM_NC 64
#endif

void
dtrlsm_simple_unblk(int m, int n, double alpha, int unitDiag,
                    const double *A, int incRowA, int incColA,
                    double *B, int incRowB, int incColB)
{
    int i, j;

    if (alpha!=1) {
        for (j=0; j<n; ++j) {
            for (i=0; i<m; ++i) {
                B[i*incRowB+j*incColB] *= alpha;
            }
        }
    }
    for (j=0; j<n; ++j) {
        dtrlsv_ref(m, unitDiag, A, incRowA, incColA, &B[j*incColB], incRowB);
    }
}

void
dtrlsm_simple_blk(int m, int n, double alpha, int unitDiag,
                  const double *A, int incRowA, int incColA,
                  double *B, int incRowB, int incColB)
{
    int i, j, l;

    const int NC = DTRLSM_NC;
    const int MC = DTRLSM_MC;

    for (j=0; j<n; j+=NC) {
        const int nc = (j+NC<=n) ? NC : n-j;

        for (i=0; i<m; i+=MC) {
            const int    mc     = (i+MC<=m) ? MC : m-i;
            const double alpha_ = (i==0) ? alpha : 1;

            dtrlsm_simple_unblk(mc, nc, alpha_, unitDiag,
                                &A[i*(incRowA+incColA)], incRowA, incColA,
                                &B[i*incRowB+j*incColB], incRowB, incColB);

            for (l=i+MC; l<m; l+= MC) {
                const int kc = (l+MC<=m) ? MC : m-l;

                dgemm(kc, nc, mc,
                      -1.0,
                      &A[l*incRowA+i*incColA], incRowA, incColA,
                      &B[i*incRowB+j*incColB], incRowB, incColB,
                      alpha_,
                      &B[l*incRowB+j*incColB], incRowB, incColB);
            }
        }
    }
}

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

#ifndef MIN_M
#define MIN_M 100
#endif

#ifndef MAX_M
#define MAX_M 1500
#endif

#ifndef INC_M
#define INC_M 100
#endif

#ifndef MIN_N
#define MIN_N 100
#endif

#ifndef MAX_N
#define MAX_N 1500
#endif

#ifndef INC_N
#define INC_N 100
#endif

#ifndef ROWMAJOR
#define ROWMAJOR 0
#endif

#ifndef ALPHA
#define ALPHA 2.5
#endif

#ifndef UNITDIAG
#define UNITDIAG 0
#endif

#if (ROWMAJOR==1)
#   define INCROW_A  MAX_M
#   define INCCOL_A  1
#else
#   define INCROW_A  1
#   define INCCOL_A  MAX_M
#endif

#if (ROWMAJOR==1)
#   define INCROW_B  MAX_N
#   define INCCOL_B  1
#else
#   define INCROW_B  1
#   define INCCOL_B  MAX_M
#endif

double A[MAX_M*MAX_M];
double B[MAX_M*MAX_M];
double X0[MAX_M*MAX_N];
double X1[MAX_M*MAX_N];

int
main()
{
    int m, n;

    randGeMatrix(MAX_M, MAX_M, A, INCROW_A, INCCOL_A);
    makeTrlDiagDom(MAX_M, UNITDIAG, A, INCROW_A, INCCOL_A);
    randGeMatrix(MAX_M, MAX_N, B, INCROW_B, INCCOL_B);

    printf("#UNITDIAG=%3d\n", UNITDIAG);
    printf("# %51s %25s\n", "dtrlsv_ref", "dtrlsv_row");
    printf("#%9s %9s %9s", "n", "INCROW_A", "INCCOL_A");
    printf(" %9s %9s", "INCROW_B", "INCCOL_B");
    printf(" %12s %12s", "t", "MFLOPS");
    printf(" %12s %12s %12s", "t", "MFLOPS", "err");
    printf("\n");

    for (m=MIN_M, n=MIN_N; m<=MAX_M && n<=MAX_N; m+=INC_M, n+=INC_N) {
        int     runs = 1;
        double  ops = (double)n*m*(m+1)/2 + (double)n*m*(m-1)/2;
        double  t0, dt, err;

        printf(" %9d %9d %9d", n, INCROW_A, INCCOL_A);
        printf(" %9d %9d", INCROW_B, INCCOL_B);

        t0   = 0;
        runs = 0;
        do {
            dgecopy(m, n, B, INCROW_B, INCCOL_B, X0, INCROW_B, INCCOL_B);
            dt = walltime();

            dtrlsm_ref(m, n, (double)ALPHA, UNITDIAG,
                       A, INCROW_A, INCCOL_A,
                       X0, INCROW_B, INCCOL_B);

            dt = walltime() - dt;
            t0 += dt;
            ++runs;
        } while (t0<0.3);
        t0 /= runs;

        printf(" %12.2e %12.2lf", t0, ops/(1000*1000*t0));

        t0   = 0;
        runs = 0;
        do {
            dgecopy(m, n, B, INCROW_B, INCCOL_B, X1, INCROW_B, INCCOL_B);
            dt = walltime();

            dtrlsm_simple_blk(m, n, (double)ALPHA, UNITDIAG,
                              A, INCROW_A, INCCOL_A,
                              X1, INCROW_B, INCCOL_B);

            dt = walltime() - dt;
            t0 += dt;
            ++runs;
        } while (t0<0.3);
        t0 /= runs;

        err = err_dtrsm(m, n, ALPHA, UNITDIAG, 1,
                        A, INCROW_A, INCCOL_A,
                        X0, INCROW_B, INCCOL_B,
                        X1, INCROW_B, INCCOL_B);

        printf(" %12.2e %12.2lf %12.2e", t0, ops/(1000*1000*t0), err);
        printf("\n");

    }

    return 0;
}