#include "libscl.h"
#include <ctime>

using namespace std;
using namespace scl;

namespace {
  bool tstequal(const realmat& a, const realmat& b);
  string dimensions(const realmat& a, const realmat& b);
  realmat crudemult(const realmat& a, const realmat& b);
  void genmat(INT_32BIT& seed, INTEGER arows, INTEGER acols, realmat& a);
}

int main()
{
  /*
  const INTEGER root_cblas_mult_size
    = INTEGER(sqrt(cblas_mult_size));
  const INTEGER fourth_root_cblas_mult_size
    = INTEGER(pow(cblas_mult_size,0.25));
  const INTEGER three_fourth_root_cblas_mult_size
    = INTEGER(pow(cblas_mult_size,0.75));
  */

  clock_t clock_start, clock_stop;

  clock_start = clock();

  INT_32BIT seed = 100542;

  cout << boolalpha;

  cout << '\t' << starbox("/Tests of gemv//");

  INTEGER arows;
  INTEGER acols;

  INTEGER brows;
  INTEGER bcols;

  realmat a,b,r1,r2;

  arows = cblas_mult_size - 5;
  acols = brows = cblas_mult_size;
  bcols = 1;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(a,b);
  dgmprd(a,b,r1);
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  
  arows = cblas_mult_size + 5;
  acols = brows = cblas_mult_size;
  bcols = 1;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(a,b);
  dgmprd(a,b,r1);
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  
  arows = cblas_mult_size - 100;
  acols = brows = 25;
  bcols = 1;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(a,b);
  dgmprd(a,b,r1);
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  
  arows = cblas_mult_size - 100;
  acols = brows = 25;
  bcols = 1;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(a,b);
  dgmprd(a,b,r1);
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  
  arows = 25;
  acols = cblas_mult_size - 10;
  brows = 1;
  bcols = 25;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(realmat(T(a)),realmat(T(b)));
  r1 = T(a)*T(b);
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  arows = 25;
  acols = cblas_mult_size + 10;
  brows = 1;
  bcols = 25;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(realmat(T(a)),realmat(T(b)));
  r1 = T(a)*T(b);
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  arows = 25;
  acols = cblas_mult_size - 10;
  brows = 25;
  bcols = 1;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(realmat(T(a)),b);
  r1 = T(a)*b;
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  arows = 25;
  acols = cblas_mult_size + 10;
  brows = 25;
  bcols = 1;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(realmat(T(a)),b);
  r1 = T(a)*b;
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  arows = cblas_mult_size - 10;
  acols = 25;
  brows = 1;
  bcols = 25;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(a,realmat(T(b)));
  r1 = a*T(b);
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  
  arows = cblas_mult_size + 10;
  acols = 25;
  brows = 1;
  bcols = 25;
  genmat(seed,arows,acols,a);
  genmat(seed,brows,bcols,b);
  r2 = crudemult(a,realmat(T(b)));
  r1 = a*T(b);
  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  clock_stop = clock();

  REAL time = REAL(clock_stop - clock_start)/REAL(CLOCKS_PER_SEC);
  cerr << "clock = " << time << '\n';

  return 0;
}

namespace {

  bool tstequal(const realmat& a, const realmat& b) {
    if (a.size() != b.size()) return false;
    for (INTEGER i=1; i<=a.size(); ++i) {
      if (fabs(a[i]-b[i]) > EPS) {
        cout << '\t' << "i = " << i << '\n';
        cout << '\t' << "a[i] = " << a[i] << '\n';
        cout << '\t' << "b[i] = " << b[i] << '\n';
        return false;
      }
    }
    return true;
  }

  string dimensions(const realmat& a, const realmat& b)
  {
    INTEGER arows = a.nrow();
    INTEGER acols = a.ncol();
    INTEGER brows = b.nrow();
    INTEGER bcols = b.ncol();
    string rv;
    rv = "arows = " + fmt('d',1,arows).get_ostr() + ", ";
    rv += "acols = " + fmt('d',1,acols).get_ostr() + ", ";
    rv += "brows = " + fmt('d',1,brows).get_ostr() + ", ";
    rv += "bcols = " + fmt('d',1,bcols).get_ostr() + ", ";
    rv += "rsize = " + fmt('d',1,arows*bcols).get_ostr() + ", ";
    return rv;
  }

  realmat crudemult(const realmat& a, const realmat& b)
  {
    INTEGER arows = a.nrow();
    INTEGER acols = a.ncol();
    INTEGER brows = b.nrow();
    INTEGER bcols = b.ncol();

    realmat r(arows,bcols,0.0);

    const REAL* ax = a.get_x();
    const REAL* bx = b.get_x();
    REAL* rx = r.get_x();

    REAL* rij = rx;
    REAL* rtop = rx + r.size();

    const REAL* aik;
    const REAL* bkj;
    for (INTEGER j=0; j<bcols; j++) {
      for (INTEGER k=0; k<acols; k++) {
        aik = ax + arows*k;
        bkj = bx + brows*j + k;
        rij = rx + arows*j;
        rtop = rij + arows;
        while (rij < rtop)
          *rij++ += *aik++ * *bkj;
      }
    }

    return r;
  }

  void genmat(INT_32BIT& seed, INTEGER arows, INTEGER acols, realmat& a)
  {
    a.resize(arows,acols);
    for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  }

}
