#include "sopnamsp.h"
#include "machdefs.h"

#include <math.h>
#include <iostream>

#include "tarrinit.h"
#include "array.h"
#include "timing.h"

/*  --------------------------------------------------------- */
/*  ---- Programme de test des operations et conversions ---  */
/*  ----    sur les TArray/TMatrix de SOPHYA    ----          */
/*  2000-2004      -   C. Magneville,  R. Ansari              */
/*  carrt -> Help/Usage                                       */
/*  Pour effectuer les tests :                                */
/* Conversions entre arr<int> et arr<float>                   */
/*  carrt  ac         -> Rc=0 OK                              */
/*  carrt  ace        -> Exception/Rc=9                       */
/*  carrt  oso        -> Rc=0 OK                              */
/*  carrt  odo        -> Rc=0 OK                              */
/*  --- carrt oso/odo NRow NCols  pour test grande tailles    */
/*  --------------------------------------------------------- */

static int test_ac(bool fgszm=false);
static int test_oso(int r, int c);
static int test_odo(int r, int c);

static int prtlev = 0;
static int nprt = 100;

template <class T>
void getMinMax(TArray<T>& a, T& min, T& max) { a.MinMax(min, max); }

int main(int narg, char* arg[])
{

  SophyaInit();
  InitTim();   // Initializing the CPU timer

  if (narg < 2) {
    cout << " carrt ac/ace/oso/odo [NRow=5] [NCols=10] [prtlev=0] [nprtmax=100] \n"
	 << "   ac : array conversion test   \n" 
	 << "   ace : array conversion test + conversion with sizemismatch -> exception/rc=9  \n" 
	 << "   oso : operation same memory organisation \n" 
	 << "   odo : operation different memory organisation \n" << endl;
    return 1;
  }

  string opt = arg[1];
  int nr = 5;
  int nc = 10;
  if (narg > 2) nr = atoi(arg[2]);
  if (narg > 3) nc = atoi(arg[3]);
  if (narg > 4) prtlev = atoi(arg[4]);
  if (narg > 5) nprt = atoi(arg[5]);

  BaseArray::SetMaxPrint(nprt);
  int rc = 0;
  PrtTim(" Start of Test ");
  try {
    if (opt == "ac") rc = test_ac(false);
    else if (opt == "ace") rc = test_ac(true);
    else if (opt == "oso") rc = test_oso(nr, nc);
    else if (opt == "odo") rc = test_odo(nr, nc);
  }
  catch (PThrowable & exc) {
    cerr << " catched Exception " << exc.Msg() << endl;
    rc = 9;
  }  
  catch (...) {
    cerr << " catched unknown (...) exception " << endl; 
  }  

  PrtTim(" End of Test ");
  cout << " ----- End of carrt - Rc= " << rc << " ----- " << endl;
  if (rc > 255) rc = 255;
  return rc;
}


int test_ac(bool fgszm)
{
  cout << "\n -----> Testing TArray Conversion <---- " << endl;
  int rc = 0;
  int_4 min, max;
  TArray<int_4> ia(7,5),iac,iad;   
  ia = RegularSequence(10., 2.);
  TArray<r_4> ra(7,5);   
  ra = ia;
  if (prtlev>0) 
    cout << " ra = ia(= RegularSequence(10., 2.)) = \n " << ra << endl;
  iac = ra;
  iad = iac-ia;
  getMinMax(iad, min, max);
  if ((min != 0) || (max != 0)) {
    cout << "(1) !!! ERROR ra(7,5)=ia(7,5)  Min=" << min << " Max=" << max << endl;
    rc += 1;
  }
  else cout << "(1) OK  ra(7,5)=ia(7,5)  OK " << endl; 

  TArray<r_4> rb(5,3);
  rb =  RegularSequence(20., .5);  
  TMatrix<int_4> mx(3,5);
  TMatrix<int_4> mxc(3,5);
  mxc = RegularSequence(20., .5);  
  if (prtlev>0) 
    cout << " TArray<r_4> rb(5,3); rb =  RegularSequence(20., .5); rb : " << rb << endl;
  mx = rb;
  if (prtlev>0) 
    cout << " TMatrix<int_4> mx(3,5);       mx = rb;  mx: " << mx << endl;
  TMatrix<int_4> mxd = mxc-mx;
  getMinMax(mxd, min, max);
  if ((min != 0) || (max != 0)) {
    cout << "(2) !!! ERROR mx<int>(5,3)=rb<r_4>(5,3)  Min=" << min << " Max=" << max << endl;
    rc += 2;
  }
  else cout << "(2) OK  mx<int>(5,3)=rb<r_4>(5,3)  OK " << endl; 

  TArray<r_4> rcv(10);
  rcv = 3.1415;
  TMatrix<r_8> mx2(rcv);
  if (prtlev>0) 
    cout << "   TArray<r_4> rcv(10) = 3.1415 , TMatrix<r_8> mx2(rc)=rcv : " << mx2 << endl;
  TMatrix<int_4> mx2d = mx2-3.1415;
  getMinMax(mxd, min, max);
  if ((min != 0) || (max != 0)) {
    cout << "(3) !!! ERROR mx2<r_8>(10)=rcv<r_4>(10)  Min=" << min << " Max=" << max << endl;
    rc += 4;
  }
  else cout << "(3) OK  mx2<r_8>(10)=rcv<r_4>(10)  OK " << endl; 
  if (fgszm) {
    TMatrix<int_4> mx3(3,5);
    cout << "(4) Trying  TMatrix<int_4> mx3(3,5);  mx3 = rcv(10); ?? -> Exception " << endl;
    mx3 = rcv;
  }
  return(rc);
}

int test_oso(int nr, int nc)
{
  cout << "\n -----> Testing TArray operation TArray<int_4>(nx=" << nr << ",ny=" 
       << nc << " ) " << endl;
  cout << " -----> OSO (Same Mem. Mapping) Same memory mapping (default) for all arrays" << endl;

  int rc = 0;
  int min,max;

  TArray<int_4> a(nr,nc);
  a = 20;
  TArray<int_4> b(nr,nc);
  b = 9;
  if (prtlev > 0) {
    cout << " A = \n " << a << endl;
    cout << " B = \n " << b << endl;
  }

  TArray<int_4> ar;
  ar = a+2;  
  if (prtlev > 0) 
    cout << "(1) AR = A+2=  " << ar << endl;
  getMinMax(ar, min, max);
  if ((min != 22) || (max != 22)) {
    cout << "(1) !!! ERROR Test AR=A+2  Min=" << min << " Max=" << max << endl;
    rc += 1;
  }
  else cout << "(1) OK   Test AR=A+2  OK " << endl; 

  ar = 28-a;  
  if (prtlev > 0) 
    cout << "(2) AR = 28-A=  " << ar << endl;
  getMinMax(ar, min, max);
  if ((min != 8) || (max != 8)) {
    cout << "(2) !!! ERROR Test AR=28-A  Min=" << min << " Max=" << max << endl;
    rc += 2;
  }
  else cout << "(2) OK   Test AR=28-A  OK " << endl;
 
  ar = a*3;  
  if (prtlev > 0) 
    cout << "(3) AR = A*3=  " << ar << endl;
  getMinMax(ar, min, max);
  if ((min != 60) || (max != 60)) {
    cout << "(3) !!! ERROR Test AR=A*3  Min=" << min << " Max=" << max << endl;
    rc += 4;
  }
  else cout << "(3) OK   Test AR=A*3  OK " << endl; 

  ar = a/5;  
  if (prtlev > 0) 
    cout << "(4) AR = A/5=  " << ar << endl;
  getMinMax(ar, min, max);
  if ((min != 4) || (max != 4)) {
    cout << "(4) !!! ERROR Test AR=A/5  Min=" << min << " Max=" << max << endl;
    rc += 8;
  }
  else cout << "(4) OK   Test AR=A/5  OK " << endl; 

  ar = 140/a;  
  if (prtlev > 0) 
    cout << "(5) AR = 140/A= \n " << ar << endl;
  getMinMax(ar, min, max);
  if ((min != 7) || (max != 7)) {
    cout << "(5) !!! ERROR Test AR=140/A  Min=" << min << " Max=" << max << endl;
    rc += 16;
  }
  else cout << "(5) OK   Test AR=140/A  OK " << endl; 

  ar -= 3;  
  if (prtlev > 0) 
    cout << "(6) AR -= 3  " << ar << endl;
  getMinMax(ar, min, max);
  if ((min != 4) || (max != 4)) {
    cout << "(6) !!! ERROR Test AR-=3  Min=" << min << " Max=" << max << endl;
    rc += 32;
  }
  else cout << "(6) OK   Test AR-=3  OK " << endl; 


  TArray<int_4> c = a-2*b;
  if (prtlev > 0) 
    cout << "(7) C = A-2*B=  " << c << endl;
  getMinMax(c, min, max);
  if ((min != 2) || (max != 2)) {
    cout << "(7) !!! ERROR Test C=A+2*B  Min=" << min << " Max=" << max << endl;
    rc += 64;
  }
  else cout << "(7) OK   Test C=A+2*B  OK " << endl; 

  (a*4+1).DivElt(b,c);
  if (prtlev > 0) cout << "(8) C = (A*4+1)/B = " << c << endl;
  getMinMax(c, min, max);
  if ((min != 9) || (max != 9)) {
    cout << "(8) !!! ERROR Test C = (A*4+1)/B  Min=" << min << " Max=" << max << endl;
    rc += 128;
  }
  else cout << "(8) OK   Test C = (A*4+1)/B ((DivElt) OK " << endl; 
  
  a = RegularSequence(1, 1);
  b = RegularSequence(nr*nc, -1);
  c = (a*2)+(b*2);
  if (prtlev > 0) {
    cout << "(9) A =  " << a << endl;
    cout << "(9) B =  " << b << endl;
    cout << "(9) C = 2A-2*B=  " << c << endl;
  }
  c /= (nr*nc+1);
  getMinMax(c, min, max);
  if ((min != 2) || (max != 2)) {
    cout << "(9) !!! ERROR Test C=2A+2*B  Min=" << min << " Max=" << max << endl;
    rc += 256;
  }
  else cout << "(9) OK   Test C=2A+2*B  OK " << endl; 

  a = 10;
  b = 4;
  c = 3*a-2*b;
  if (prtlev > 0)  cout << "(10) A=10, B=4, C = 3*A-2*B = " << c << endl;
  getMinMax(c, min, max);
  if ((min != 22) || (max != 22)) {
    cout << "(10) !!! ERROR Test C = 3*A-2*B  Min=" << min << " Max=" << max << endl;
    rc += 512;
  }
  else cout << "(10) OK   Test C = 3*A-2*B  OK " << endl; 

  c = 3*a+2*b;
  if (prtlev > 0)  cout << "(11) A=10, B=4, C = 3*A+2*B = " << c << endl;
  getMinMax(c, min, max);
  if ((min != 38) || (max != 38)) {
    cout << "(11) !!! ERROR Test C = 3*A+2*B  Min=" << min << " Max=" << max << endl;
    rc += 512;
  }
  else cout << "(11) OK   Test C = 3*A+2*B  OK " << endl; 

  a.MulElt(b+2,c);
  if (prtlev > 0)  cout << "(12) A=10, B=4, a.MulElt(b+2,c) = " << c << endl;
  getMinMax(c, min, max);
  if ((min != 60) || (max != 60)) {
    cout << "(12) !!! ERROR Test a.MulElt(b+2,c) Min=" << min << " Max=" << max << endl;
    rc += 1024;
  }
  else cout << "(12) OK   Test A.MulElt(B+2,C)  OK " << endl; 

  (3*a).DivElt(b+2,c);
  if (prtlev > 0)  cout << "(13) A=10, B=4, (3*a).DivElt(b+2,c) C = " << c << endl;
  getMinMax(c, min, max);
  if ((min != 5) || (max != 5)) {
    cout << "(13) !!! ERROR Test (3*a).DivElt(b+2,c) Min=" << min << " Max=" << max << endl;
    rc += 1024;
  }
  else cout << "(13) OK  Test (3*A).DivElt(B+2,C)  OK " << endl; 
  
  a = 3;
  b = 2;
  int_4 sc = a.ScalarProduct(b);
  if (prtlev > 0)  cout << "(14) A=3, B=3, a.ScalarProduct(b)= " << sc << endl;
  if (sc != nr*nc*6)   {
    cout << "(14) !!! ERROR Test  a.ScalarProduct(b) sc " << sc << " !=" << nr*nc*6 << endl;
    rc += 2048;
  }
  else cout << "(14) OK  Test a.ScalarProduct(b)  OK " << endl; 

  TMatrix<int_4> mx(nr, nc), mxc(nr,nc);
  for(int ir=0; ir<nr; ir++)
    for(int ic=0; ic<nc; ic++) mxc(ir, ic) = ir*1000+ic*6;
  if (prtlev > 0)  cout << "(15) mx.Row()/Column() operations MXC=" << mxc << endl;
  
  TVector<int_4> cv(nr), rv(nc, BaseArray::RowVector);
  cv = RegularSequence(0, 1);
  rv = RegularSequence(0, 1);
  for(int ir=0; ir<nr; ir++)  { 
    mx.Row(ir) = rv*3;
    mx.Row(ir) *= 2;
  }
  for(int ic=0; ic<nc; ic++)  { 
    mx.Column(ic) += cv*1000;
  }
  TMatrix<int_4> mxd = mx-mxc;
  if (prtlev > 0)  { 
    cout << "(15)  MX=" << mx << endl;
    cout << "(15)  MXD=" << mxd << endl;
  }
  getMinMax(mxd, min, max);
  if ((min != 0) || (max != 0)) {
    cout << "(15) !!! ERROR Test mx.Row()/Column() operations Min=" << min << " Max=" << max << endl;
    rc += 8192;
  }
  else cout << "(15) OK  Test mx.Row()/Column() operations  OK " << endl; 

  a = RegularSequence(1, 1);
  b = -a; 
  c = a+b;
  if (prtlev > 0) {
    cout << "(16) A =  " << a << endl;
    cout << "(16) B =  ((? == -A)" << b << endl;
    cout << "(16) C = A+B  " << c << endl;
  }
  getMinMax(c, min, max);
  if ((min != 0) || (max != 0)) {
    cout << "(16) !!! ERROR Test B=-A  Min=" << min << " Max=" << max << endl;
    rc += 256;
  }
  else cout << "(16) OK   Test B=-A  OK " << endl; 
  
  return(rc);
}

int test_odo(int nr, int nc)
{
  cout << "\n -----> Testing TMatrix operation TMatrix<r_4>(nr=" << nr << ",nc=" 
       << nc << " ) " << endl;
  cout << " --> ODO (Diff. Mem. Mapping) A,CC: CMemoryMapping AR,B,C: FortranMemoryMapping" << endl;

  int rc = 0;
  r_4 min,max;

  TMatrix<r_4> a(nr,nc,BaseArray::CMemoryMapping);
  a = 20;
  TMatrix<r_4> b(nr,nc,BaseArray::FortranMemoryMapping);
  b = 9;

  TMatrix<r_4> ar(nr,nc,BaseArray::FortranMemoryMapping);

  if (prtlev > 0) {
    cout << " A = \n " << a << endl;
    cout << " B = \n " << b << endl;
  }

  a.AddCst(2.,ar);
  if (prtlev > 0) 
    cout << "(1) AR = A+2 : a.AddCst(2,ar) : " << ar << endl;
  getMinMax(ar, min, max);
  if ((fabs(min-22.) > 0.0001) || (fabs(max-22.) > 0.0001)) {
    cout << "(1) !!! ERROR Test AR=A+2  a.AddCst(2,ar)   Min=" << min << " Max=" << max << endl;
    rc += 1;
  }
  else cout << "(1) OK   Test AR=A+2 a.AddCst(2,ar)  OK " << endl; 

  a.SubCst(28.,ar,true);
  if (prtlev > 0) 
    cout << "(2) AR = 28-A : a.SubCst(28,ar,true) :  " << ar << endl;
  getMinMax(ar, min, max);
  if ((fabs(min-8.) > 0.0001) || (fabs(max-8.) > 0.0001)) {
    cout << "(2) !!! ERROR Test AR=28-A  a.SubCst(28,ar,true) Min=" << min << " Max=" << max << endl;
    rc += 2;
  }
  else cout << "(2) OK   Test AR=28-A  a.SubCst(28,ar,true) OK " << endl;
 
  a.MulCst(3.,ar);
  if (prtlev > 0) 
    cout << "(3) AR = A*3= a.MulCst(3,ar)  " << ar << endl;
  getMinMax(ar, min, max);
  if ((fabs(min-60.) > 0.0001) || (fabs(max-60.) > 0.0001)) {
    cout << "(3) !!! ERROR Test AR=A*3  a.MulCst(3,ar) Min=" << min << " Max=" << max << endl;
    rc += 4;
  }
  else cout << "(3) OK   Test AR=A*3  a.MulCst(3,ar) OK " << endl; 

  a.DivCst(5.,ar);
  if (prtlev > 0) 
    cout << "(4) AR = A/5= a.DivCst(5,ar) : " << ar << endl;
  getMinMax(ar, min, max);
  if ((fabs(min-4.) > 0.0001) || (fabs(max-4.) > 0.0001)) {
    cout << "(4) !!! ERROR Test AR=A/5  a.DivCst(5,ar) Min=" << min << " Max=" << max << endl;
    rc += 8;
  }
  else cout << "(4) OK   Test AR=A/5  a.DivCst(5,ar) OK " << endl; 

  a.DivCst(140.,ar,true);
  if (prtlev > 0) 
    cout << "(5) AR = 140/A= a.DivCst(140.,ar,true) : " << ar << endl;
  getMinMax(ar, min, max);
  if ((fabs(min-7.) > 0.0001) || (fabs(max-7.) > 0.0001)) {
    cout << "(5) !!! ERROR Test AR=140/A  a.DivCst(140.,ar,true) Min=" << min << " Max=" << max << endl;
    rc += 16;
  }
  else cout << "(5) OK   Test AR=140/A  a.DivCst(140.,ar,true) OK " << endl; 

  ar -= 3.;  
  if (prtlev > 0) 
    cout << "(6) AR -= 3 " << ar << endl;
  getMinMax(ar, min, max);
  if ((min != 4) || (max != 4)) {
    cout << "(6) !!! ERROR Test AR-=3  Min=" << min << " Max=" << max << endl;
    rc += 32;
  }
  else cout << "(6) OK   Test AR-=3  OK " << endl; 

  TMatrix<r_4> c(nr,nc,BaseArray::FortranMemoryMapping);
  c = a-(b*2.0f);
  if (prtlev > 0) cout << "(7) C = A-2*B=  " << c << endl;

  getMinMax(c, min, max);
  if ((fabs(min-2.) > 0.0001) || (fabs(max-2.) > 0.0001)) {
    cout << "(7) !!! ERROR Test C=A+2*B  Min=" << min << " Max=" << max << endl;
    rc += 64;
  }
  else cout << "(7) OK   Test C=A+2*B  OK " << endl; 

  (a*4.0f+5.5f).DivElt(b,c);
  if (prtlev > 0) cout << "(8) C = (A*4+12)/B = " << c << endl;
  getMinMax(c, min, max);
  if ((fabs(min-9.5) > 0.0001) || (fabs(max-9.5) > 0.0001)) {
    cout << "(8) !!! ERROR Test C = (A*4+12)/B  Min=" << min << " Max=" << max << endl;
    rc += 128;
  }
  else cout << "(8) OK   Test C = (A*4+12)/B  OK " << endl; 

  a = RegularSequence(1, 1);
  b = RegularSequence(nr*nc, -1);
  c = (a*2.f)+(b*2.f);
  if (prtlev > 0) {
    cout << "(9) A =  " << a << endl;
    cout << "(9) B =  " << b << endl;
    cout << "(9) C = 2A-2*B= " << c << endl;
  }
  c /= (nr*nc+1);
  getMinMax(c, min, max);
  if ((fabs(min-2.) > 0.0001) || (fabs(max-2.) > 0.0001)) {
    cout << "(9) !!! ERROR Test C=2A+2*B  Min=" << min << " Max=" << max << endl;
    rc += 256;
  }
  else cout << "(9) OK   Test3 C=2A+2*B  OK " << endl; 

// ----------------------
// Operation de type array +-*/ array 

  a = 10.;
  b = 4.;

  (3.f*a).SubElt(2.f*b,c);
  if (prtlev > 0)  cout << "(10) A=10, B=4, (3*a).SubElt(2*b,c)  C=" << c << endl;
  getMinMax(c, min, max);
  if ((fabs(min-22.) > 0.0001) || (fabs(max-22.) > 0.0001)) {
    cout << "(10) !!! ERROR Test (3*a).SubElt(2*b,c)  Min=" << min << " Max=" << max << endl;
    rc += 512;
  }
  else cout << "(10) OK   Test (3*a).SubElt(2.*b,c)  OK " << endl; 

  (3.f*a).AddElt(2.f*b,c);
  if (prtlev > 0)  cout << "(11) A=10, B=4, (3*a).AddElt(2*b,c) C= " << c << endl;
  getMinMax(c, min, max);
  if ((fabs(min-38.) > 0.0001) || (fabs(max-38.) > 0.0001)) {
    cout << "(11) !!! ERROR Test (3*A).AddElt(2.*B,C)  Min=" << min << " Max=" << max << endl;
    rc += 512;
  }
  else cout << "(11) OK   Test (3*a).AddElt(2*b,c) OK " << endl; 

  a.MulElt(b+2.f,c);
  if (prtlev > 0)  cout << "(12) A=10, B=4, a.MulElt(b+2,c) = " << c << endl;
  getMinMax(c, min, max);
  if ((fabs(min-60.) > 0.0001) || (fabs(max-60.) > 0.0001)) {
    cout << "(12) !!! ERROR Test a.MulElt(b+2,c) Min=" << min << " Max=" << max << endl;
    rc += 1024;
  }
  else cout << "(12) OK   Test A.MulElt(B+2,C)  OK " << endl; 

  (3.f*a).DivElt(b+2.f,c);
  if (prtlev > 0)  cout << "(13) A=10, B=4, (3*a).DivElt(b+2,c) C = " << c << endl;
  getMinMax(c, min, max);
  if ((fabs(min-5.) > 0.0001) || (fabs(max-5.) > 0.0001)) {
    cout << "(13) !!! ERROR Test (3*a).DivElt(b+2,c) Min=" << min << " Max=" << max << endl;
    rc += 1024;
  }
  else cout << "(13) OK  Test (3*A).DivElt(B+2,C)  OK " << endl; 

  TMatrix<r_4> cc(nr,nc,BaseArray::CMemoryMapping);

  (3.f*a).SubElt(2.f*b,cc);
  if (prtlev > 0)  cout << "(14) A=10, B=4, (3*a).SubElt(2*b,cc)  " << cc << endl;
  getMinMax(cc, min, max);
  if ((fabs(min-22.) > 0.0001) || (fabs(max-22.) > 0.0001)) {
    cout << "(14) !!! ERROR Test (3*a).SubElt(2*b,cc)  Min=" << min << " Max=" << max << endl;
    rc += 2048;
  }
  else cout << "(14) OK   Test (3*A).SubElt(2*B,CC)  OK " << endl; 


  (3.f*a).AddElt(2.f*b,cc);
  if (prtlev > 0)  cout << "(15) A=10, B=4, (3*a).AddElt(2*b,cc) CC= " << cc << endl;
  getMinMax(cc, min, max);
  if ((fabs(min-38.) > 0.0001) || (fabs(max-38.) > 0.0001)) {
    cout << "(15) !!! ERROR Test (3*a).AddElt(2*b,c)  Min=" << min << " Max=" << max << endl;
    rc += 4096;
  }
  else cout << "(15) OK   Test (3*A).AddElt(2*B,CC) OK " << endl; 

  a.MulElt(b+2.f,cc);
  if (prtlev > 0)  cout << "(16) A=10, B=4, a.MulElt(b+2,c) = " << cc << endl;
  getMinMax(cc, min, max);
  if ((fabs(min-60.) > 0.0001) || (fabs(max-60.) > 0.0001)) {
    cout << "(16) !!! ERROR Test A.MulElt(B+2,CC) Min=" << min << " Max=" << max << endl;
    rc += 8192;
  }
  else cout << "(16) OK   Test A.MulElt(B+2,C)  OK " << endl; 

  (3.f*a).DivElt(b+2.f,cc);
  if (prtlev > 0)  cout << "(17) A=10, B=4, (3*a).DivElt(b+2,cc) C = " << cc << endl;
  getMinMax(cc, min, max);
  if ((fabs(min-5.) > 0.0001) || (fabs(max-5.) > 0.0001)) {
    cout << "(17) !!! ERROR Test (3*a).DivElt(b+2,c) Min=" << min << " Max=" << max << endl;
    rc += 16384;
  }
  else cout << "(17) OK  Test (3*A).DivElt(B+2,CC)  OK " << endl; 

  a = 3.1;
  b = -2.2;
  r_4 sc = a.ScalarProduct(b);
  r_4 sce = 0.;
  for(int kkk=0; kkk<nr*nc; kkk++) sce += -2.2*3.1;
  if (prtlev > 0)  cout << "(18) A=3, B=3, a.ScalarProduct(b)= " << sc << endl;
  sc -= sce;
  if ((fabs(sc) > 0.00001) || (fabs(sc) > 0.00001)) {
    cout << "(18) !!! ERROR Test  a.ScalarProduct(b) sc-" <<  sce
	 << " = " << sc << endl;
    rc += 32768;
  }
  else cout << "(18) OK  Test a.ScalarProduct(b)  OK " << endl; 

  TMatrix<r_4> mx(nr, nc), mxc(nr,nc,BaseArray::FortranMemoryMapping);
  for(int ir=0; ir<nr; ir++)
    for(int ic=0; ic<nc; ic++) mxc(ir, ic) = ir*1000+ic*6;
  if (prtlev > 0)  cout << "(19) mx.Row()/Column() operations MXC=" << mxc << endl;
  
  TVector<r_4> cv(nr), rv(nc, BaseArray::RowVector, BaseArray::FortranMemoryMapping);
  cv = RegularSequence(0, 1);
  rv = RegularSequence(0, 1);
  for(int ir=0; ir<nr; ir++)  { 
    mx.Row(ir) = rv*3.f;
    mx.Row(ir) *= 2.f;
  }
  for(int ic=0; ic<nc; ic++)  { 
    mx.Column(ic) += cv*1000.f;
  }
  TMatrix<r_4> mxd = mx-mxc;
  if (prtlev > 0)  { 
    cout << "(19)  MX=" << mx << endl;
    cout << "(19)  MXD=" << mxd << endl;
  }
  getMinMax(mxd, min, max);
  if ((fabs(min) > 0.0001) || (fabs(max) > 0.0001)) {
    cout << "(19) !!! ERROR Test mx.Row()/Column() operations Min=" << min << " Max=" << max << endl;
    rc += 65536;
  }
  else cout << "(19) OK  Test mx.Row()/Column() operations  OK " << endl; 

  a = RegularSequence(3.14, 1.3);
  b = -a; 
  c = a+b;
  if (prtlev > 0) {
    cout << "(20) A =  " << a << endl;
    cout << "(20) B =  ((? == -A)" << b << endl;
    cout << "(16) C = A+B  " << c << endl;
  }
  getMinMax(c, min, max);
  if ((fabs(min) > 0.0001) || (fabs(max) > 0.0001)) {
    cout << "(20) !!! ERROR Test B=-A  Min=" << min << " Max=" << max << endl;
    rc += 131072;
  }
  else cout << "(20) OK   Test B=-A  OK " << endl; 
  
  return(rc);

  return(rc);

}