#include "scltypes.h"
#include "realmat.h"
#include "sclfuncs.h"
#include <stdio.h>
#include <ctime>

using namespace std;
using namespace scl;

namespace {
  bool tstequal(const realmat& a, const realmat& b);
  realmat tstindx(const realmat& a, realmat& b);
  void loop(const realmat& a, const realmat& b, realmat& r);
}

LIB_ERROR_HANDLER_PTR previous; 

void do_not_exit (string msg)
{
  throw lib_error(msg);  
}

int main()
{

  clock_t clock_start, clock_stop;
  clock_start = clock();

  ofstream of03("ft03f001.fba");
  if ( !of03 ) {
    cerr << "Error, main, cannot open ft03f001.fba";
    exit(1);
  } 
  ostream & os03=of03;
  
  INTEGER i,j;
  INTEGER rr, cc;

  rr=6;
  cc=5;

  realmat a;

  a.resize(rr,cc,1);

  for(i=1; i<=6; i++)
  for(j=1; j<=5; j++)
    a(i,j)=i+j;

  realmat b = a;
  cout << starbox("/check copy constructor//");
  if (a == b) cout << "Copy constructor works\n";

  b = tstindx(a,b);
  cout << starbox("/does [] work?///");
  if (a == b) cout << "It does work\n";

  cout << starbox("/a followed by T(a) then by realmat ww = T(a) + T(a)//");
  //cout << a << T(a); 
  realmat ww = T(a) + realmat(T(a));
  cout << ww;

  cout << starbox("/a has elements i+j//");
  cout << a;
  os03 << starbox("/a has elements i+j//");
  os03 << a;

  cout << starbox("/Check self assignment/aa=a,a=a,a==aa/aa+=aa+=aa/bb=aa=a//");
  realmat aa=a;
  a=a;
  cout << "\n     (a == aa) = " << (a == aa) << '\n';
  aa+=aa+=aa;
  cout << aa;
  realmat bb;
  bb=aa=a;
  cout << bb;

  cout << starbox("\nCheck division\na/5.0==a/5\n(1.0/5.0)*a\na/5.0\n\n",'\n');
  cout << "\n     (a/5.0 == a/5) = " << (a/REAL(5.0) == a/5) <<'\n';
  cout << (1.0/5.0)*a << a/REAL(5.0) ;

  cout << starbox("\nCheck division\nT(a)/5.0==T(a)/5\n(1.0/5.0)*T(a)\nT(a)/5.0\n\n",'\n');
  cout << "\n     (T(a)/5.0==T(a)/5) = " << (T(a)/REAL(5.0)==T(a)/5) <<'\n';
  cout << (1.0/5.0)*T(a) << T(a)/REAL(5.0);

  cout << starbox("/T(a)*a//");
  cout << T(a)*a;
  cout << realmat(T(a))*a;
  cout << starbox("/invpsd(T(a)*a)//");
  cout << invpsd(T(a)*a);
  cout << starbox("/inv(T(a)*a)//");
  cout << inv(T(a)*a,1.e-10);
  cout << starbox("/a*invpsd(T(a)*a)*T(a)*a//");
  cout << a*invpsd(T(a)*a)*T(a)*a;
  cout << starbox("/a*inv(T(a)*a)*T(a)*a//");
  cout << a*inv(T(a)*a)*T(a)*a;
  cout << starbox("/a*ginv(T(a)*a)*T(a)*a//");
  cout << a*ginv(T(a)*a)*T(a)*a;

  INTEGER rank;
  cout << starbox("/T(a)*a//");
  cout << T(a)*a;
  cout << starbox("/inv(T(a)*a,rank)//");
  cout << inv(T(a)*a,rank);
  cout << "\n     rank = " << rank << '\n';
  cout << starbox("/ginv(T(a)*a,rank)//");
  cout << ginv(T(a)*a,rank);
  cout << "\n     rank = " << rank << '\n';
  cout << starbox("/invpsd(T(a)*a,rank)//");
  cout << invpsd(T(a)*a,rank);
  cout << "\n     rank = " << rank << '\n';

  const REAL five = 5.0;

  cout << starbox("/five*T(a)*a//");
  cout << five*T(a)*a;
  cout << starbox("/T(a)*a*five//");
  cout << T(a)*a*five;
  cout << starbox("/a*invpsd(five*T(a)*a)*T(a)*a*five//");
  cout << a*invpsd(five*T(a)*a)*T(a)*a*five;
  cout << starbox("/a*inv(T(a)*a*five)*five*T(a)*a//");
  cout << a*inv(T(a)*a*five)*five*T(a)*a;

  cout << starbox("/Inner Product//");
  realmat left = a("1"," ");
  realmat rite = realmat(T(left));
  cout << "\n     left " << left;
  cout << "\n     rite " << rite;
  cout << "\n     left*rite"    << left*rite;
  cout << "\n     left*T(left)" << left*T(left);
  cout << "\n     T(rite)*rite" << T(rite)*rite;

  cout << starbox("/a transpose, three ways:/cout<<T(a)/T(a)(i,j)/T(a)[i]//");
  cout << T(a);
  cout << '\n';
  for(i=1; i<=5; i++) {
  for(j=1; j<=6; j++) {
    cout << ' ' << T(a)(i,j); 
    if (j==6) cout << '\n';
  }
  }
  cout << '\n';
  for (i=1; i<=a.size(); i++) {
    cout << ' ' << T(a)[i] << '\n'; 
    if (!(i%5)) cout << '\n';
  }

  cout << starbox("/Ta=T(a)//");
  realmat Ta;
  Ta = T(a);
  cout << Ta;
  cout << starbox("/T(a)//");
  cout << T(a);
  cout << starbox("/Ta*a//");
  cout << Ta*a;
  cout << starbox("/T(a)*a//");
  cout << T(a)*a;
  cout << starbox("/a*Ta//");
  cout << a*Ta;
  cout << starbox("/a*T(a)//");
  cout << a*T(a);
  cout << a*realmat(T(a));

  for(j=1; j<=5; j++)
  for(i=1; i<=6; i++)
    if (i<=j) 
      a(i,j)=1; 
    else 
      a(i,j)=0;
  cout << starbox("/a = upper triangle of ones//");
  cout << a;

  cout << starbox("/T(a)*a//");
  cout << T(a)*a;
  cout << starbox("/invpsd(T(a)*a)//");
  cout << invpsd(T(a)*a);
  cout << starbox("/T(a)*a*invpsd(T(a)*a)//");
  cout << T(a)*a*invpsd(T(a)*a);

  b.resize(6,5);
  b=2.0*a;
  cout << starbox("/b=2.0*a//");
  cout << b;

  realmat c = a+b;
  cout << starbox("/c = a+b//");
  cout << c;
  cout << "\n\tc has " << c.nrow() <<" rows and " << c.ncol() << " columns.\n";

  cout << starbox("/T(c)//");
  cout << T(c);
  cout << "\n\tT(c) has " << T(c).nrow() << " rows and " 
       << T(c).ncol() << " columns.\n";
  cout << "\n\tT(c) has size " << T(c).size() << "\n.";
  cout << "\n\tT(c) has " << T(c).get_rows() << " rows and " 
       << T(c).get_cols() << " columns.\n";
  cout << "\n\tT(c) has size " << T(c).size() << "\n.";

  cout << starbox("/c=c-c//");
  c = c-c;
  cout << c;

  realmat d(6,6,2);
  cout << starbox("/d(6,6,2)/d*a//");
  cout << d*a;

  realmat e(6,6,-1);
  cout << starbox("/e(6,6,-1)/e*a//");
  cout << e*a;

  os03.flush();

  cout << starbox("/-e and +e//");
  cout << -e << +e;

  cout << starbox("/a += a//");  
  a += a;
  cout << a;
  cout << starbox("/a -= a//"); 
  a -= a;
  cout << a;

  cout << starbox("/++a  then --a//"); 
//cout << ++a << --a;  Both MS CL and Sun g++ get this wrong.
  ++a; cout << a; --a; cout << a;  //This works with them all.
  
  cout << starbox("//line 1/line 2/line 3///");
  cout << starbox("\n\nline 1\nline 2\nline 3\n\n\n",'\n');

  realmat z;
  realmat zz(5,5,0.0);
  realmat zzz(5,5);
  cout << starbox("/zz=z where z is null/before and after =//");
  cout << zz;
  zz = z;
  cout << zz;

  previous = set_lib_error_handler(&do_not_exit);
  cout << starbox("/Try to invert a null matrix//");
  try {
    invpsd(zz);
  }  
  catch (lib_error x) { 
    cout << x.msg;
  }
  previous = set_lib_error_handler(previous);
    
  previous = set_lib_error_handler(&do_not_exit);
  cout << starbox("/Try to multiply two null's//");
  try {
    zzz = zz * z;
    cout << zzz;
  }
  catch (lib_error x) { 
    cout << x.msg;
  }
  previous = set_lib_error_handler(previous);
 
  previous = set_lib_error_handler(&do_not_exit);
  cout << starbox("/Try to multiply nonconformable matrices//");
  z.resize(3,3);
  zz.resize(5,5);
  try {
    zzz = zz * z;
  }  
  catch (lib_error x) { 
    cout << x.msg;
  }
  previous = set_lib_error_handler(previous);

  cout << starbox("/test unary operators/e = -e; e = (-2)*e; e = (0.5)*e; e = +e; /e before and after//");

  cout << e;

  e = -e; e = (-2)*e; e = (0.5)*e; e = +e;

  cout << e;

  realmat nullmat;

  cout << starbox("/e and a//") << e << a;

  cout << starbox("/e == e, e != e, a == e, a != e, nullmat == nullmat//");

  cout <<" "<< (e==e) <<" "<< (e!=e) <<" "<< (a==e) 
       <<" "<< (a!=e) <<" "<< (nullmat==nullmat);

  cout << starbox("/e <= e, e < e, a <= e, a < e,/nullmat <= nullmat, nullmat < nullmat//");

  cout <<" "<< (e<=e) <<" "<< (e<e) <<" "<< (a<=e) <<" "<< (a<e)
       <<" "<< (nullmat<=nullmat) <<" "<< (nullmat<nullmat);

  cout << starbox("/e >= e, e > e, a >= e, a > e,/nullmat >= nullmat, nullmat > nullmat//");

  cout <<" "<< (e>=e) <<" "<< (e>e) <<" "<< (a>=e) <<" "<< (a>e)
       <<" "<< (nullmat>=nullmat) <<" "<< (nullmat>nullmat);
 
  cout << starbox("/e <= -e, e < -e, e >= -e, e > -e//");

  cout <<" "<< (e <= -e) <<" "<< (e < -e) <<" "<< (e >= -e)  <<" "<< (e > -e);

  previous = set_lib_error_handler(&do_not_exit);
  cout << starbox("/Try to increment a null matrix//");
  try {
    ++nullmat;
  }  
  catch (lib_error x) { 
    cout << x.msg;
  }
  previous = set_lib_error_handler(previous);

  previous = set_lib_error_handler(&do_not_exit);
  cout << starbox("/Try to decrement a null matrix//");
  try {
    --nullmat;
  }  
  catch (lib_error x) { 
    cout << x.msg;
  }
  previous = set_lib_error_handler(previous);

  cout << starbox("/Check precedence ++a - a then a++ - a//"); 
  z = a;  zz = ++z - a;  cout << zz << z; 
  z = a;  zz = z++ - a;  cout << zz << z;  cout << '\n';

  cout << starbox("/Check fill and get_x//");
  fill(a); cout << a;
  fill(a,5.0); dgmpnt(cout, a.get_x(), a.nrow(), a.ncol()) << a;

  cout << starbox("/The matrix d and submatrix s(sr,sc)//");
  for (i=1; i<=d.nrow()*d.ncol(); i++)
    d[i]=i;
  const char* sr = "1,2,4:6"; 
  const char* sc = "1:3,5:6,1,1,1";
  cout << d <<'\n'<< "     sr = " << sr << "  sc = " << sc << d(sr,sc) <<'\n';

  previous = set_lib_error_handler(&do_not_exit);
  cout << starbox("/Try a bad string/The matrix d and submatrix s(sr,sc)//");
  for (i=1; i<=d.nrow()*d.ncol(); i++)
    d[i]=i/5;
  d=d*5;
  const char* sr0 = "1,2,4:6"; 
  const char* sc0 = "1:3,5:6-,8";
  //char* sc0 = "1:3,5:2147483647,8";
  //char* sc0 = "1:3,5: -2147483648,8";
  cout << d <<'\n'<< "     sr0 = " << sr0 << "  sc0 = " << sc0 <<'\n';
  realmat w;  
  try {
    w = d(sr0,sc0);
  }
  catch (lib_error x) { 
    cout << x.msg;
  }
  previous = set_lib_error_handler(previous);
  cout << w <<'\n';

  for (i=1; i<=d.nrow()*d.ncol(); i++)  d[i]=i;
  cout << starbox("/d//") << d;

  const char* s0 = "1,3,4 ";
  const char* s1 = "1:2   ";
  const char* s2 = "2:1   "; 
  const char* s4 = "2:4   ";
  const char* s5 = "5     ";

  cout << "s0, s1, s2, s4, s5 = " << s0 << s1 << s2 << s4 << s5 <<'\n';

  cout << starbox("/cbind(d(s0,s1),d(s0,s2))//") 
                 << cbind(d(s0,s1),d(s0,s2));

  cout << starbox("/rbind(d(s4,s0),d(s5,s0))//") 
                 << rbind(d(s4,s0),d(s5,s0));

  cout << starbox("/How's this?//");
  realmat a11(3,3,1.0); 
  realmat a21(3,3,0.0);
  realmat a12=a21;
  realmat a22=a11;
  a = rbind(cbind(a11,a12),cbind(a21,a22));
  cout << a << a("1:2,5:6","1:2,5:6");

  cout << starbox("/Check subscripting//");
  for (i=1; i<=d.size(); i++) {
    d.check(i)=i;
    cout << d[i] << ' ' << d.check(i) << ' ' << d.check(i) << endl;
    cout << T(d)[i] << ' ' << T(d).check(i) << endl; 
  }
  for (i=1; i<=d.nrow(); i++) {
    for (j=1; j<=d.ncol(); j++) {
      d.check(i,j)=d.nrow()*(j-1)+i;
      cout << d(i,j) << ' ' << d.check(i,j) << ' ' << d.check(i,j) << endl;
      cout << T(d)(i,j) << ' ' << T(d).check(i,j) << endl;
    }
  }

  cout << starbox("/Here is d once again, then jvec and ivec//") << d << endl;
  
  intvec jvec = seq(2,4);
  intvec ivec = seq(-5,7);

  for (INTEGER j=1; j<=jvec.size(); j++) cout << "  " << jvec[j];
  cout << endl;
  for (INTEGER i=1; i<=ivec.size(); i++) cout << "  " << ivec[i];
  cout << endl;

  cout << starbox("/Here is d(ivec,jvec) //") << d(ivec,jvec);

  cout << starbox("/Here is d once again, then jvec and ivec//") << d << endl;
  
  jvec = seq(-2,-4);
  ivec = seq(-5,7);

  for (INTEGER j=1; j<=jvec.size(); j++) cout << "  " << jvec[j];
  cout << endl;
  for (INTEGER i=1; i<=ivec.size(); i++) cout << "  " << ivec[i];
  cout << endl;

  cout << starbox("/Here is d(ivec,jvec) //") << d(ivec,jvec);

  cout << starbox("/Here is d once again, then jvec and ivec//") << d << endl;
  
  jvec.resize(4);

  jvec[1] = 2; 
  jvec[2] = 2; 
  jvec[3] = 1; 
  jvec[4] = 1; 

  ivec = seq(-5,7);

  for (INTEGER j=1; j<=jvec.size(); j++) cout << "  " << jvec[j];
  cout << endl;
  for (INTEGER i=1; i<=ivec.size(); i++) cout << "  " << ivec[i];
  cout << endl;

  cout << starbox("/Here is d(ivec,jvec) //") << d(ivec,jvec);

  cout << starbox("/First d, then d's rows and columns//") << d;
  for (INTEGER i=1; i<=d.nrow(); i++) cout << d(i,"");
  for (INTEGER j=1; j<=d.ncol(); j++) cout << d("",j);

  starbox("/a = T(a); a = T(a);//");
  for (INTEGER i=1; i<=a.size(); ++i) a[i] = i;
  cout << a; 
  a = T(a);
  cout << a; 
  a = T(a);
  cout << a - a;

  cout << starbox("/check a*realmat(T(b))//");
  b = a("1,2","");
  //cout << a << b;
  cout << a*T(b);
  cout << a*realmat(T(b));
  
  cout << starbox("/check realmat(T(a))*b//");
  b = a("","1,2");
  //cout << a << b;
  cout << T(a)*b;
  cout << realmat(T(a))*b;
  
  cout << starbox("/check realmat_cmp//");
  realmat_cmp compare;
  realmat lhs(1,15,1.0);
  realmat rhs(1,15,1.0);
  cout << boolalpha << compare(lhs,rhs) << '\n';
  rhs++;
  cout << compare(lhs,rhs) << '\n';
  rhs--; rhs--;
  cout << compare(lhs,rhs) << '\n';

  cout << starbox("/check resize to null matrix//");
  a.resize(0,0);
  cout << a;

  cout << starbox("/check push_back()//");
  realmat p;
  for (INTEGER i=1; i<=5; ++i) p.push_back(REAL(i));
  cout << p;
  p.resize(0,0);
  for (INTEGER i=1; i<=6; ++i) p.push_back(REAL(i));
  cout << T(p);
  for (INTEGER i=1; i<=6; ++i) p.push_back(REAL(i));
  cout << T(p);

  rr=6;
  cc=5;

  a.resize(rr,cc);
  b.resize(cc,rr);

  for(INTEGER j=1; j<=cc; j++) {
    for(INTEGER i=1; i<=rr; i++) {
      a(i,j)=i+j;
      b(j,i)=i+j;
    }
  }

  cout << starbox("/check T(a) + T(a)//");
  //cout << a << b;
  cout << T(a) + T(a) << realmat(T(a)) + realmat(T(a)) 
       << T(a) + T(a) -  (realmat(T(a)) + realmat(T(a))) ;
        
  cout << starbox("/check T(a) + b//");
  //cout << a << b;
  cout << T(a) + b  << realmat(T(a)) + b 
       << T(a) + b -  (realmat(T(a)) + b) ;
        
  cout << starbox("/check a + T(b)//");
  //cout << a << b;
  cout << a + T(b)  << realmat(T(a)) + b 
       << a + T(b) -  (a + realmat(T(b))) ;


  cout << starbox("/check  +T(b)//");
  //cout << a << b;
  cout << +T(b)  << +realmat(T(b)) 
       << +T(b)  - +realmat(T(b));

  cout << starbox("/check T(a) - T(a)//");
  //cout << a << b;
  cout << T(a) - T(a) << realmat(T(a)) - realmat(T(a)) 
       << T(a) - T(a) -  (realmat(T(a)) - realmat(T(a))) ;
        
  cout << starbox("/check T(a) - b//");
  //cout << a << b;
  cout << T(a) - b  << realmat(T(a)) - b 
       << T(a) - b -  (realmat(T(a)) - b) ;
        
  cout << starbox("/check a - T(b)//");
  //cout << a << b;
  cout << a - T(b)  << a - realmat(T(b))
       << a - T(b) -  (a - realmat(T(b))) ;


  cout << starbox("/check  -T(b)//");
  //cout << a << b;
  cout << -T(b)  << -realmat(T(b)) 
       << -T(b)  - (-realmat(T(b)));

        
  cout << starbox("/check T(a)*T(b)//");
  //cout << a << b;
  cout << T(a)*T(b) << T(b*a) << T(a)*T(b) - T(b*a);

  cout << starbox("/check T(a)*b//");
  b = T(b);
  //cout << a << b;
  cout << T(a)*b << realmat(T(a))*b << T(a)*b - realmat(T(a))*b;

  cout << starbox("/check a*T(b)//");
  //cout << a << b;
  cout << a*T(b) << a*realmat(T(b)) << a*T(b) - a*realmat(T(b));

  cout << starbox("/check iterators//");
  a.resize(25,30,8.0);
  b.resize(25,30);
  realmat::iterator aij = a.begin();
  realmat::const_iterator bij = b.begin();
  while (aij != a.end()) *aij++ = *bij++;
  if (a == b) cout << "\n\t iterators work" << '\n';

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

  realmat r0;
  realmat r1;

  cout << '\t' << '1' << '\n';
  a.resize(cblas_mult_size+5,1,0.25);
  b.resize(1,cblas_mult_size+5,0.25);
  loop(a,b,r0);
  r1 = a*b;
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  realmat c0 = r0;
  cout << '\t' << tstequal(r0,c0) << '\n'; 
  c0 = r0;
  cout << '\t' << tstequal(r0,c0) << '\n'; 

  cout << '\t' << '2' << '\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);
  loop(a,b,r0);
  r1 = a*b;
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  realmat c1 = r0;
  cout << '\t' << tstequal(r0,c1) << '\n'; 
  c1 = r0;
  cout << '\t' << tstequal(r0,c1) << '\n'; 

  cout << '\t' << '3' << '\n';
  a.resize(5,1,0.25);
  b.resize(1,5,0.25);
  loop(a,b,r0);
  r1 = a*b;
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  realmat c2 = r0;
  cout << '\t' << tstequal(r0,c2) << '\n'; 
  c2 = r0;
  cout << '\t' << tstequal(r0,c2) << '\n'; 

  cout << '\t' << '4' << '\n';
  a.resize(5,5,0.25);
  b.resize(5,5,0.25);
  loop(a,b,r0);
  r1 = a*b;
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  realmat c3 = r0;
  cout << '\t' << tstequal(r0,c3) << '\n'; 
  c3 = r0;
  cout << '\t' << tstequal(r0,c3) << '\n'; 

  cout << '\t' << '5' << '\n';
  a.resize(1,1,0.25);
  b.resize(1,1,0.25);
  loop(a,b,r0);
  r1 = a*b;
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  realmat c4 = r0;
  cout << '\t' << tstequal(r0,c4) << '\n'; 
  c4 = r0;
  cout << '\t' << tstequal(r0,c4) << '\n'; 

  cout << '\t' << '6' << '\n';
  a.resize(1,1,0.25);
  b.resize(1,1,0.25);
  loop(a,b,r0);
  r1 = a*b;
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  realmat c5 = r0;
  cout << '\t' << tstequal(r0,c5) << '\n'; 
  c5 = r0;
  cout << '\t' << tstequal(r0,c5) << '\n'; 

  INT_32BIT seed = 5;

  cout << '\t' << '7' << '\n';
  a.resize(cblas_mult_size+15,cblas_mult_size+5);
  b.resize(cblas_mult_size+5,cblas_mult_size+15);
  for (INTEGER j=1; j<=a.ncol(); ++j) {
    for (INTEGER i=1; i<=a.nrow(); ++i) {
      a.check(i,j) = unsk(seed);
      b.check(j,i) = a(i,j);
    }
  }
  loop(b,a,r0);
  r1 = T(a)*a;
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  loop(a,b,r0);
  r1 = a*T(a);
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  loop(b,a,r0);
  r1 = T(a)*T(b);
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  c = T(a);
  r1 = c*a;
  cout << '\t' << tstequal(r0,r1) << '\n'; 

  cout << '\t' << '8' << '\n';
  a.resize(15,5);
  b.resize(5,15);
  for (INTEGER j=1; j<=a.ncol(); ++j) {
    for (INTEGER i=1; i<=a.nrow(); ++i) {
      a.check(i,j) = unsk(seed);
      b.check(j,i) = a(i,j);
    }
  }
  loop(b,a,r0);
  r1 = T(a)*a;
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  loop(a,b,r0);
  r1 = a*T(a);
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  loop(b,a,r0);
  r1 = T(a)*T(b);
  cout << '\t' << tstequal(r0,r1) << '\n'; 
  c = T(a);
  r1 = c*a;
  cout << '\t' << tstequal(r0,r1) << '\n'; 

  clock_stop = clock();
  std::cerr << "clock = "
    << REAL(clock_stop - clock_start)/REAL(CLOCKS_PER_SEC) << '\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]) > fabs(a[i])*1.0e-10) {
        cout << '\t' << "i = " << i << '\n';
        cout << '\t' << "a[i] = " << fmt('f',26,17,a[i]) << '\n';
        cout << '\t' << "b[i] = " << fmt('f',26,17,b[i]) << '\n';
        return false;
      }
    }
    return true;
  }

  realmat tstindx(const realmat& a, realmat& b)
  {
    b.resize(a.nrow(),a.ncol());
    for (INTEGER i=1; i<= a.size(); ++i) {
      b[i] = a[i];
    }
    return b;
  }

  void loop(const realmat& a, const realmat& b, realmat& r) 
  {
    INTEGER arows = a.get_rows();
    INTEGER acols = a.get_cols();
    INTEGER brows = b.get_rows();
    INTEGER bcols = b.get_cols();
  
    if (acols != brows) error("Error, loop, matrices not conformable."); 
    if (a.size() == 0) error("Error, loop, null matrix.");
  
    const REAL zero = 0.0;
  
    INTEGER newrows = arows;  
    INTEGER newcols = bcols;
    INTEGER newlen = newrows*newcols;
  
    if (r.get_rows() != newrows || r.get_cols() != newcols) {
      r.resize(newrows,newcols);
    }
  
    REAL* newx = r.get_x();
    const REAL* ax = a.get_x();
    const REAL* bx = b.get_x();
  
    if (newlen == 1) {   //Special case of an inner product.
  
      const REAL* ai = ax;
      const REAL* bi = bx;
      const REAL* btop = bi + brows;  
      *newx = zero;
      while (bi < btop) *newx += *ai++ * *bi++;
  
    } 
    else {               //General case.
  
      REAL* rij = newx;
      REAL* rtop = newx + newlen;
      while (rij < rtop) *rij++ = zero;
  
      const REAL* aik;
      const REAL* bkj;
      for (INTEGER j=0; j<newcols; j++) {
        for (INTEGER k=0; k<acols; k++) {
          aik = ax + arows*k;
          bkj = bx + brows*j + k;
          rij = newx + newrows*j;
          rtop = rij + newrows;
          while (rij < rtop)
            *rij++ += *aik++ * *bkj;
        }
      }
    }
  }

}