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

using namespace std;
using namespace scl;

namespace {
  bool tstequal(const realmat& a, const realmat& b);
  string dimensions(realmat& a, realmat& b);
}

int main()
{

  const INTEGER loop_count = 500;

  cout << boolalpha;

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

  INTEGER arows;
  INTEGER acols;

  INTEGER brows;
  INTEGER bcols;

  realmat a,b,r1,r2;

  arows = brows = acols = bcols = 1;

  a.resize(arows,acols);
  b.resize(brows,bcols);

  INT_32BIT seed = 100542;

  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);

  dgmprd(a,b,r1);

  r2 = a*b;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,b,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  
  //dgmprd(realmat(), b, r1);
  //dgmprd(a, realmat(), r1);


  arows = 1 ; acols = 25 ; brows = 25; bcols = 1;

  a.resize(arows,acols);
  b.resize(brows,bcols);

  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);

  dgmprd(a,b,r1);

  r2 = a*b;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,b,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';
  
  arows = 5 ; acols = 25 ; brows = 25; bcols = 3;

  a.resize(arows,acols);
  b.resize(brows,bcols);

  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);

  dgmprd(a,b,r1);

  r2 = a*b;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,b,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  intvec aiv, ajv, biv, bjv;

  aiv = seq(1,arows);
  ajv = seq(1,acols);
  
  biv = seq(1,brows);
  bjv = seq(1,bcols);

  realmat aa, bb;

  dgmprd(a,aiv,ajv,b,biv,bjv,r1);

  r2 = a(aiv,biv)*b(biv,bjv);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,aiv,ajv,b,biv,bjv,r1);

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(aa,b,biv,bjv,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,aiv,ajv,bb,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  aiv[2] = 0;
  bjv[2] = 0;

  dgmprd(a,aiv,ajv,b,biv,bjv,r1);

  r2 = a(aiv,ajv)*b(biv,bjv);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,aiv,ajv,b,biv,bjv,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(aa,b,biv,bjv,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,aiv,ajv,bb,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';


  ajv.resize(1,2);
  biv.resize(1,1);

  dgmprd(a,aiv,ajv,b,biv,bjv,r1);

  r2 = a(aiv,ajv)*b(biv,bjv);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,aiv,ajv,b,biv,bjv,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(aa,b,biv,bjv,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,aiv,ajv,bb,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';


  aiv = seq(1,arows);
  ajv = seq(1,acols);
  
  biv = seq(1,brows);
  bjv = seq(1,bcols);

  aiv.resize(1,1);
  ajv[acols] = -1;

  biv[brows] = -1;
  bjv.resize(1,2);

  dgmprd(a,aiv,ajv,b,biv,bjv,r1);

  r2 = a(aiv,ajv)*b(biv,bjv);

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,aiv,ajv,b,biv,bjv,r1);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(aa,bb,r1);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(aa,b,biv,bjv,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  dgmprd(a,aiv,ajv,bb,r1);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  cout << '\t' << starbox("/cblas_mult_sizes crossover timing//");

  clock_t clock_start, clock_stop;

  INTEGER lim = INTEGER(sqrt(cblas_mult_size));

  arows = lim+1 ; acols = lim+1 ; brows = lim+1; bcols = lim+1;

  a.resize(arows,acols);
  b.resize(brows,bcols);

  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);

  clock_start = clock();

  aiv = seq(1,arows);
  ajv = seq(1,acols);
  
  biv = seq(1,brows);
  bjv = seq(1,bcols);

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,b,biv,bjv,r1);
  }

  r2 = a(aiv,biv)*b(biv,bjv);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,b,biv,bjv,r1);
  }

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(aa,b,biv,bjv,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,bb,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  aiv[2] = 0;
  bjv[2] = 0;

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,b,biv,bjv,r1);
  }

  r2 = a(aiv,ajv)*b(biv,bjv);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,b,biv,bjv,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(aa,b,biv,bjv,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,bb,r1);
  }

  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';

  arows = lim-1 ; acols = lim-1 ; brows = lim-1; bcols = lim-1;

  a.resize(arows,acols);
  b.resize(brows,bcols);

  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);

  clock_start = clock();

  aiv = seq(1,arows);
  ajv = seq(1,acols);
  
  biv = seq(1,brows);
  bjv = seq(1,bcols);

  for (INTEGER i=1; i<=loop_count; ++i) {
   dgmprd(a,aiv,ajv,b,biv,bjv,r1);
  }

  r2 = a(aiv,biv)*b(biv,bjv);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,b,biv,bjv,r1);
  }

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(aa,b,biv,bjv,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,bb,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  aiv[2] = 0;
  bjv[2] = 0;

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,b,biv,bjv,r1);
  }

  r2 = a(aiv,ajv)*b(biv,bjv);

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,b,biv,bjv,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  aa = a(aiv,ajv);

  bb = b(biv,bjv);

  r2 = aa*bb;

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(aa,b,biv,bjv,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,aiv,ajv,bb,r1);
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  clock_stop = clock();

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

  const INTEGER tst = 100; 

  arows = tst+3 ; acols = tst+5 ; brows = tst+5; bcols = tst+3;

  a.resize(arows,acols);
  b.resize(brows,bcols);

  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);
  
  clock_start = clock();

  for (INTEGER i=1; i<=loop_count; ++i) {
    dgmprd(a,b,r1);
    r2 = a*b;
  }

  cout << '\t' << dimensions(a,b) << tstequal(r1,r2) << '\n';

  clock_stop = clock();

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

  cout << '\t' << starbox("/checks on matrices over and under cblas sizes//");

  realmat r0;

  a.resize(cblas_mult_size+5,1,0.25);
  b.resize(1,cblas_mult_size+5,0.25);
  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);
  dgmprd(a,b,r0);
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r0,r1) << '\n'; 
  realmat c0 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c0) << '\n'; 
  c0 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c0) << '\n'; 

  a.resize(cblas_mult_size+5,cblas_mult_size+5,0.25);
  b.resize(cblas_mult_size+5,cblas_mult_size+5,0.25);
  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);
  dgmprd(a,b,r0);
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r0,r1) << '\n'; 
  realmat c1 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c1) << '\n'; 
  c1 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c1) << '\n'; 

  a.resize(5,1,0.25);
  b.resize(1,5,0.25);
  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);
  dgmprd(a,b,r0);
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r0,r1) << '\n'; 
  realmat c2 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c2) << '\n'; 
  c2 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c2) << '\n'; 

  a.resize(5,5,0.25);
  b.resize(5,5,0.25);
  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);
  dgmprd(a,b,r0);
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r0,r1) << '\n'; 
  realmat c3 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c3) << '\n'; 
  c3 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c3) << '\n'; 

  a.resize(1,1,0.25);
  b.resize(1,1,0.25);
  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);
  dgmprd(a,b,r0);
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r0,r1) << '\n'; 
  realmat c4 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c4) << '\n'; 
  c4 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c4) << '\n'; 

  a.resize(1,1,0.25);
  b.resize(1,1,0.25);
  for (INTEGER i=1; i<=a.size(); ++i) a[i] = unsk(seed);
  for (INTEGER i=1; i<=b.size(); ++i) b[i] = unsk(seed);
  dgmprd(a,b,r0);
  r1 = a*b;
  cout << '\t' << dimensions(a,b) << tstequal(r0,r1) << '\n'; 
  realmat c5 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c5) << '\n'; 
  c5 = r0;
  cout << '\t' << dimensions(a,b) << tstequal(r0,c5) << '\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(realmat& a, 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() + ", ";
    return rv;
  }
}
