/* --- Programme de test des matrices carres speciales
   --- This code is part of the SOPHYA library ---
  (C) Univ. Paris-Sud   (C) LAL-IN2P3/CNRS   (C) IRFU-CEA
  (C) R. Ansari, C.Magneville    2009  

   Appel: 
   # Pour verification visuelle (/impression)
   csh> tssqmx 0     
   # Test matrices diagonales (Rc=0 -> OK) 
   csh> tssqmx d 
   # Test matrices triangulaires (Rc=0 -> OK) 
   csh> tssqmx t 
   # Test matrices symmetriques (Rc=0 -> OK) 
   csh> tssqmx s 
---------------------------------------------------------- */



#include <iostream>
#include <string>

#include "trngmtx.h"
#include "diagmtx.h"
#include "symmtx.h"
#include "fiospsqmtx.h"
#include "tarrinit.h"

using namespace std;
using namespace SOPHYA;

// --- Declaration des fonctions de ce fichier 
static int SSQM_Check0();
static int TriangMtxCheck(sa_size_t NR=7);
static int DiagMtxCheck(sa_size_t NR=7);
static int SymmMtxCheck(sa_size_t NR=7);


static int PrtLev = 0;  
//-----------------------------------------------------------
//-----------------------------------------------------------
int main(int narg, char* arg[])
{
int rc = 0;
if(narg<2) {
  cout << "Usage: tssqmx Action \n Action= 0/t/d/s [NR=7] [PrtLev=0] " << endl;
  return 1;
}
sa_size_t NR = 7;
if (narg>2)  NR = atoi(arg[2]);
if (narg>3)  PrtLev = atoi(arg[3]);
try {
// appel enregistrement des handlers PPF (devra etre fait par l'initialiseur de TArray
  SophyaInit();
  cout << " ---- tssqmx : Test of Special Square Matrix classes ----- " << endl;
  if (*arg[1] == 'd')  rc = DiagMtxCheck(NR);
  else if (*arg[1] == 't')  rc = TriangMtxCheck(NR);
  else if (*arg[1] == 's')  rc = SymmMtxCheck(NR);
  else  rc = SSQM_Check0();
}
catch (PThrowable exc) {
  cerr << " tssqmx.cc: catched Exception " << exc.Msg() << endl;
  rc = 97;
}  
catch (...) {
  cerr << " tssqmx.cc: catched unknown (...) exception " << endl; 
  rc  = 99;
}  

cout << " ---- End of tssqmx : Rc= " << rc << " ---- " << endl;
return rc;
}



//-------------------------
//  Fonction de test avec impression 
static int SSQM_Check0()
{
  cout << " ========== SSQM_Check0() =========== " << endl;
  cout << "1- LowerTriangularMatrix<int_4> trmx(5) ..." << endl;
  LowerTriangularMatrix<int_4> trmx(5);
  cout << "2-  trmx = 5 ..." << endl;
  //  int_4 c = 5;
  trmx = 5;
  //  trmx.SetCst(6);
  cout << "3-  trmx += 3 ..." << endl;
  trmx += 3;
  cout << "4- cout << trmx , t2=trmx t2.print " << endl; 
  cout << trmx; 
  LowerTriangularMatrix<int_4> t2;
  t2 = trmx;
  t2.Print(cout, 10);

  cout << "5.a- tt = trmx + t2 " << endl; 
  LowerTriangularMatrix<int_4> tt;
  tt = trmx + t2;
  cout << tt ; 
  tt.Print(cout, 60);
  cout << "5.b- tt2 = tt - 4*t2 " << endl; 
  LowerTriangularMatrix<int_4> tt2 = tt - 4*t2;
  tt2.Print(cout, 60);

  cout << "6- LowerTriangularMatrix<r_4> rtt; rtt= RegularSequence(5.6, 0.2)" << endl;
  LowerTriangularMatrix<r_4> rtt(7);
  rtt= RegularSequence(5.6, 0.2);
  rtt.Print(cout,60);
  
  cout << "11- diag = 16 " << endl;
  DiagonalMatrix<int_4> diag(5), dmx(7);
  diag = 16;
  cout << diag; 
  diag.Print(cout);

  cout << "11.b- dmx =  RegularSequence(1., 2.) " << endl;
  dmx = RegularSequence(1., 2.);
  dmx.Print(cout);
  cout << "12- dd = diag*2+3 , dd2 = dd/5-diag/4+4 " << endl;
  DiagonalMatrix<int_4> dd = diag*2+3;
  dd.Print(cout);
  DiagonalMatrix<int_4> dd2;
  dd2 = dd/5-diag/4+4;
  dd2.Print(cout);

  {
  cout << "12- Writing dd , dd2 to POutPersist po(diagmtx.ppf) " << endl;
  POutPersist po("diagmtx.ppf");
  po << dd << dd2;
  }
  {
  cout << "12.b- Reading rdd , rdd2 from PInPersist po(diagmtx.ppf) " << endl;
  PInPersist pi("diagmtx.ppf");
  DiagonalMatrix<int_4> rdd,rdd2; 
  pi >> rdd >> rdd2;
  rdd.Print(cout);
  rdd2.Print(cout);
  }


  cout << "21- Symmetric = RegularSequence " << endl;
  SymmetricMatrix<int_4> sym(5);
  sym = RegularSequence(1,2);
  sym.Print(cout,-1);

  {
  cout << "22- Writing sym to POutPersist po(symmtx.ppf) " << endl;
  POutPersist po("symmtx.ppf");
  po << sym ;
  }
  {
  cout << "22.b- Reading rsym from PInPersist po(symmtx.ppf) " << endl;
  PInPersist pi("symmtx.ppf");
  SymmetricMatrix<int_4> rsym; 
  pi >> rsym;
  rsym.Print(cout,-1);
  }

  cout << "23: Matrix multiplication on SymmetricMatrix<r_8>" << endl;
  SymmetricMatrix<r_8>  tmxa(3), tmxb(3);
  tmxa = RandomSequence(RandomSequence::Gaussian);
  tmxb = RandomSequence(RandomSequence::Gaussian,3.);
  TMatrix<r_8> mxd = tmxa * tmxb;
  tmxa.Print(cout,-1);
  tmxb.Print(cout,-1);
  mxd.Print(cout,-1);

  return 0;  
}


//-------------------------
//  Fonction de test pour matrices triangulaires
static int TriangMtxCheck(sa_size_t NR)
{
  cout << " ========== TriangMtxCheck() LowerTriangularMatrix<T> check - NR= " << NR 
       << " =========== " << endl;
  int rc = 0;
//----- Test 1
  cout << "1/ Test 1 : Element Access Check LowerTriangularMatrix<int_8>" << endl;
  LowerTriangularMatrix<int_8>  tmx(NR), tmx2(NR);
  tmx2 = RegularSequence(5., 3.);
  int_8 vv=5;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=c; r<tmx.NRows(); r++) {
      tmx(r,c) = vv; vv+=3;
    } 	  
  }
  int nzz = 0;
  for (sa_size_t r=0; r<tmx.NRows()-1; r++) {
    for (sa_size_t c=r+1; c<tmx.NCols(); c++) {
      tmx(r,c) = 4.5; if (tmx(r,c)!=0) nzz++;
    } 	  
  }
  if (nzz>0) {
    cout << " !!! Test_1.a ERROR : nzz= " << nzz << endl; 
    rc += 1;
  }
  else cout << " ---> Test_1.a OK " << endl;

  LowerTriangularMatrix<int_8> td8 = tmx-tmx2;
  int_8 min8, max8;
  td8.MinMax(min8, max8);
  if (PrtLev>2) {
    cout << " LowerTriangularMatrix<int_8> tmx : " << endl;
    tmx.Print(cout, -1); 
    cout << " LowerTriangularMatrix<int_8> td8=tmx-tmx2 : " << endl;
    td8.Print(cout, -1); 
  }
  if ((min8!=0)||(max8!=0)) {
    cout << " !!! Test_1.b ERROR : diff::min=" << min8 << " diff::max=" << max8 << endl;
    rc += 2;
  }
  else cout << " ---> Test_1.b OK " << endl;

//----- Test 2
  cout << "2/ Test 2 : Operator =, PPF write check LowerTriangularMatrix<int_8>" << endl;
  {
  LowerTriangularMatrix<int_8> ctmx;
  ctmx = tmx;
  POutPersist po("ttrngmtx.ppf");
  po << ctmx;
  }
  {
  PInPersist pi("ttrngmtx.ppf");
  LowerTriangularMatrix<int_8> rtmx; 
  pi >> rtmx;
  td8 = rtmx-tmx;
  if (PrtLev>2) {
    cout << " LowerTriangularMatrix<int_8> rtmx : " << endl;
    rtmx.Print(cout, -1); 
  }
  }
  td8.MinMax(min8, max8);
  if (PrtLev>2) {
    cout << " LowerTriangularMatrix<int_8> td8=tmx-tmx2 : " << endl;
    td8.Print(cout, -1); 
  }
  if ((min8!=0)||(max8!=0)) {
    cout << " !!! Test_2 ERROR : diff::min=" << min8 << " diff::max=" << max8 << endl;
    rc += 4;
  }
  else cout << " ---> Test_2 OK " << endl;


//----- Test 3
  cout << "3/ Test 3 : ConvertToStdmatrix() + operations on LowerTriangularMatrix<r_8>" << endl;
  LowerTriangularMatrix<r_8>  tdmx(NR);
  tdmx = RandomSequence(RandomSequence::Gaussian);
  TMatrix<r_8> dmx = tdmx.ConvertToStdMatrix();
  int nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (tdmx(r,c) != dmx(r,c))  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_3.a ERROR : nerr = " << nerr << endl;
    rc += 8;
  }
  else cout << " ---> Test_3.a OK " << endl;


  tdmx *= M_PI;    tdmx -= 0.67;
  dmx *= M_PI;    dmx -= 0.67;
  LowerTriangularMatrix<r_8>  tdmx2(dmx);
  LowerTriangularMatrix<r_8> tdd8 = tdmx-tdmx2;
  r_8 mind, maxd;
  tdd8.MinMax(mind, maxd);
  if (PrtLev>2) {
    cout << " LowerTriangularMatrix<r_8> tdmx : " << endl;
    tdmx.Print(cout, -1); 
    cout << " LowerTriangularMatrix<int_8> tdd8=tdmx-tdmx2 : " << endl;
    tdd8.Print(cout, -1); 
  }
  if ((fabs(mind)>1.e-19)||(fabs(maxd)>1e-19)) {
    cout << " !!! Test_3.b ERROR : diff::min=" << mind << " diff::max=" << maxd << endl;
    rc += 16;
  }
  else cout << " ---> Test_3.b OK " << endl;

  LowerTriangularMatrix<r_8>  tdmxa(NR);
  tdmxa = RegularSequence(2.,3.);
  TMatrix<r_8> dmxa = tdmxa.ConvertToStdMatrix();
  LowerTriangularMatrix<r_8>  tdmxb = (74.32*((5.89*tdmx-tdmxa)&&tdmxa)+14.5)/tdmx;
  TMatrix<r_8>  dmxb = (74.32*((5.89*dmx-dmxa)&&dmxa)+14.5)/dmx;
  LowerTriangularMatrix<r_8>  tdmxc(dmxb);
  tdd8 = tdmxb-tdmxc;
  tdd8.MinMax(mind, maxd);
  if (PrtLev>2) {
    cout << " LowerTriangularMatrix<r_8> tdmxb : " << endl;
    tdmx.Print(cout, -1); 
    cout << " LowerTriangularMatrix<int_8> tdd8=tdmxb-tdmxc : " << endl;
    tdd8.Print(cout, -1); 
  }
  if ((fabs(mind)>1.e-19)||(fabs(maxd)>1e-19)) {
    cout << " !!! Test_3.c ERROR : diff::min=" << mind << " diff::max=" << maxd << endl;
    rc += 32;
  }
  else cout << " ---> Test_3.c OK " << endl;

//----- Test 4
  cout << "4/ Test 4 : Matrix multiplication on LowerTriangularMatrix<r_8>" << endl;
  LowerTriangularMatrix<r_8>  tmxa(NR), tmxb(NR);
  tmxa = RandomSequence(RandomSequence::Gaussian);
  tmxb = RandomSequence(RandomSequence::Gaussian,3.);
  LowerTriangularMatrix<r_8>  tmxc = tmxa * tmxb;
  TMatrix<r_8> mxa = tmxa.ConvertToStdMatrix();
  TMatrix<r_8> mxb = tmxb.ConvertToStdMatrix();
  TMatrix<r_8> mxc = mxa * mxb;

  nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (fabs(tmxc(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.a ERROR : nerr = " << nerr << endl;
    rc += 64;
  }
  else cout << " ---> Test_4.a OK " << endl;

  mxb = RandomSequence(RandomSequence::Gaussian, 9.);
  TMatrix<r_8> mxd = tmxa*mxb;
  mxc = mxa*mxb;
  if (PrtLev>2) {
    cout << " DiagonalMatrix<r_8> tmxa : " << endl;
    tmxa.Print(cout, -1); 
    cout << " TMatrix<r_8> mxb : " << endl;
    mxb.Print(cout, -1); 
    cout << " TMatrix<r_8> mxd=tmxa*mxb : " << endl;
    mxd.Print(cout, -1); 
  }

  nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (fabs(mxd(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.b ERROR : nerr = " << nerr << endl;
    rc += 128;
  }
  else cout << " ---> Test_4.b OK " << endl;

  mxd = mxb*tmxa;
  mxc = mxb*mxa;
  nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (fabs(mxd(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.c ERROR : nerr = " << nerr << endl;
    rc += 256;
  }
  else cout << " ---> Test_4.c OK " << endl;


  return rc;
}

//-------------------------
//  Fonction de test pour matrices diagonales
static int DiagMtxCheck(sa_size_t NR)
{
  cout << " ========== DiagMtxCheck() DiagonalMatrix<T> check - NR= " << NR 
       << " =========== " << endl;
  int rc = 0;
//----- Test 1
  cout << "1/ Test 1 : Element Access Check DiagonalMatrix<int_8>" << endl;
  DiagonalMatrix<int_8>  tmx(NR), tmx2(NR);
  tmx2 = RegularSequence(5., 3.);
  int_8 vv=5;
  for (sa_size_t r=0; r<tmx.NRows(); r++) {
    tmx(r,r) = vv; vv+=3;
  }
  int nzz = 0;
  for (sa_size_t r=0; r<tmx.NRows(); r++) {
    for (sa_size_t c=0; c<tmx.NCols(); c++) {
      if (r==c) continue;
      tmx(r,c) = 4.5; if (tmx(r,c)!=0) nzz++;
    } 	  
  }
  if (nzz>0) {
    cout << " !!! Test_1.a ERROR : nzz= " << nzz << endl; 
    rc += 1;
  }
  else cout << " ---> Test_1.a OK " << endl;

  DiagonalMatrix<int_8> td8 = tmx-tmx2;
  int_8 min8, max8;
  td8.MinMax(min8, max8);
  if (PrtLev>2) {
    cout << " DiagonalMatrix<int_8> tmx : " << endl;
    tmx.Print(cout, -1); 
    cout << " DiagonalMatrix<int_8> td8=tmx-tmx2 : " << endl;
    td8.Print(cout, -1); 
  }
  if ((min8!=0)||(max8!=0)) {
    cout << " !!! Test_1.b ERROR : diff::min=" << min8 << " diff::max=" << max8 << endl;
    rc += 2;
  }
  else cout << " ---> Test_1.b OK " << endl;

//----- Test 2
  cout << "2/ Test 2 : Operator =, PPF write check DiagonalMatrix<int_8>" << endl;
  {
  DiagonalMatrix<int_8> ctmx;
  ctmx = tmx;
  POutPersist po("ttrngmtx.ppf");
  po << ctmx;
  }
  {
  PInPersist pi("ttrngmtx.ppf");
  DiagonalMatrix<int_8> rtmx; 
  pi >> rtmx;
  td8 = rtmx-tmx;
  if (PrtLev>2) {
    cout << " DiagonalMatrix<int_8> rtmx : " << endl;
    rtmx.Print(cout, -1); 
  }
  }
  td8.MinMax(min8, max8);
  if (PrtLev>2) {
    cout << " DiagonalMatrix<int_8> td8=tmx-tmx2 : " << endl;
    td8.Print(cout, -1); 
  }
  if ((min8!=0)||(max8!=0)) {
    cout << " !!! Test_2 ERROR : diff::min=" << min8 << " diff::max=" << max8 << endl;
    rc += 4;
  }
  else cout << " ---> Test_2 OK " << endl;


//----- Test 3
  cout << "3/ Test 3 : ConvertToStdmatrix() + operations on DiagonalMatrix<r_8>" << endl;
  DiagonalMatrix<r_8>  tdmx(NR);
  tdmx = RandomSequence(RandomSequence::Gaussian);
  TMatrix<r_8> dmx = tdmx.ConvertToStdMatrix();
  int nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (tdmx(r,c) != dmx(r,c))  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_3.a ERROR : nerr = " << nerr << endl;
    rc += 8;
  }
  else cout << " ---> Test_3.a OK " << endl;


  tdmx *= M_PI;    tdmx -= 0.67;
  dmx *= M_PI;    dmx -= 0.67;
  DiagonalMatrix<r_8>  tdmx2(dmx);
  DiagonalMatrix<r_8> tdd8 = tdmx-tdmx2;
  r_8 mind, maxd;
  tdd8.MinMax(mind, maxd);
  if (PrtLev>2) {
    cout << " DiagonalMatrix<r_8> tdmx : " << endl;
    tdmx.Print(cout, -1); 
    cout << " DiagonalMatrix<int_8> tdd8=tdmx-tdmx2 : " << endl;
    tdd8.Print(cout, -1); 
  }
  if ((fabs(mind)>1.e-19)||(fabs(maxd)>1e-19)) {
    cout << " !!! Test_3.b ERROR : diff::min=" << mind << " diff::max=" << maxd << endl;
    rc += 16;
  }
  else cout << " ---> Test_3.b OK " << endl;

  DiagonalMatrix<r_8>  tdmxa(NR);
  tdmxa = RegularSequence(2.,3.);
  TMatrix<r_8> dmxa = tdmxa.ConvertToStdMatrix();
  DiagonalMatrix<r_8>  tdmxb = (74.32*((5.89*tdmx-tdmxa)&&tdmxa)+14.5)/tdmx;
  TMatrix<r_8>  dmxb = (74.32*((5.89*dmx-dmxa)&&dmxa)+14.5)/dmx;
  DiagonalMatrix<r_8>  tdmxc(dmxb);
  tdd8 = tdmxb-tdmxc;
  tdd8.MinMax(mind, maxd);
  if (PrtLev>2) {
    cout << " DiagonalMatrix<r_8> tdmxb : " << endl;
    tdmx.Print(cout, -1); 
    cout << " DiagonalMatrix<int_8> tdd8=tdmxb-tdmxc : " << endl;
    tdd8.Print(cout, -1); 
  }
  if ((fabs(mind)>1.e-19)||(fabs(maxd)>1e-19)) {
    cout << " !!! Test_3.c ERROR : diff::min=" << mind << " diff::max=" << maxd << endl;
    rc += 32;
  }
  else cout << " ---> Test_3.c OK " << endl;

//----- Test 4
  cout << "4/ Test 4 : Matrix multiplication on DiagonalMatrix<r_8>" << endl;
  DiagonalMatrix<r_8>  tmxa(NR), tmxb(NR);
  tmxa = RandomSequence(RandomSequence::Gaussian);
  tmxb = RandomSequence(RandomSequence::Gaussian,3.);
  DiagonalMatrix<r_8>  tmxc = tmxa * tmxb;
  TMatrix<r_8> mxa = tmxa.ConvertToStdMatrix();
  TMatrix<r_8> mxb = tmxb.ConvertToStdMatrix();
  TMatrix<r_8> mxc = mxa * mxb;

  nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (fabs(tmxc(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.a ERROR : nerr = " << nerr << endl;
    rc += 64;
  }
  else cout << " ---> Test_4.a OK " << endl;

  mxb = RandomSequence(RandomSequence::Gaussian, 9.);
  TMatrix<r_8> mxd = tmxa*mxb;
  mxc = mxa*mxb;
  if (PrtLev>2) {
    cout << " DiagonalMatrix<r_8> tmxa : " << endl;
    tmxa.Print(cout, -1); 
    cout << " TMatrix<r_8> mxb : " << endl;
    mxb.Print(cout, -1); 
    cout << " TMatrix<r_8> mxd=tmxa*mxb : " << endl;
    mxd.Print(cout, -1); 
  }

  nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (fabs(mxd(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.b ERROR : nerr = " << nerr << endl;
    rc += 128;
  }
  else cout << " ---> Test_4.b OK " << endl;

  mxd = mxb*tmxa;
  mxc = mxb*mxa;
  nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (fabs(mxd(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.c ERROR : nerr = " << nerr << endl;
    rc += 256;
  }
  else cout << " ---> Test_4.c OK " << endl;


  return rc;
}

//-------------------------
//  Fonction de test pour matrices diagonales
static int SymmMtxCheck(sa_size_t NR)
{
  cout << " ========== SymmMtxCheck() SymmetricMatrix<T> check - NR= " << NR 
       << " =========== " << endl;
  int rc = 0;
//----- Test 1
  cout << "1/ Test 1 : Element Access Check SymmetricMatrix<int_8>" << endl;
  SymmetricMatrix<int_8>  tmx(NR), tmx2(NR);
  tmx2 = RegularSequence(5., 3.);
  int_8 vv=5;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=c; r<tmx.NRows(); r++) {
      tmx(r,c) = vv; vv+=3;
    } 	  
  }
  int nzz = 0;
  for (sa_size_t r=0; r<tmx.NRows()-1; r++) {
    for (sa_size_t c=r+1; c<tmx.NCols(); c++) {
      if (tmx(r,c)!=tmx(c,r)) nzz++;
    } 	  
  }
  if (nzz>0) {
    cout << " !!! Test_1.a ERROR : nzz= " << nzz << endl; 
    rc += 1;
  }
  else cout << " ---> Test_1.a OK " << endl;

  SymmetricMatrix<int_8> td8 = tmx-tmx2;
  int_8 min8, max8;
  td8.MinMax(min8, max8);
  if (PrtLev>2) {
    cout << " SymmetricMatrix<int_8> tmx : " << endl;
    tmx.Print(cout, -1); 
    cout << " SymmetricMatrix<int_8> td8=tmx-tmx2 : " << endl;
    td8.Print(cout, -1); 
  }
  if ((min8!=0)||(max8!=0)) {
    cout << " !!! Test_1.b ERROR : diff::min=" << min8 << " diff::max=" << max8 << endl;
    rc += 2;
  }
  else cout << " ---> Test_1.b OK " << endl;

//----- Test 2
  cout << "2/ Test 2 : Operator =, PPF write check SymmetricMatrix<int_8>" << endl;
  {
  SymmetricMatrix<int_8> ctmx;
  ctmx = tmx;
  POutPersist po("ttrngmtx.ppf");
  po << ctmx;
  }
  {
  PInPersist pi("ttrngmtx.ppf");
  SymmetricMatrix<int_8> rtmx; 
  pi >> rtmx;
  td8 = rtmx-tmx;
  if (PrtLev>2) {
    cout << " SymmetricMatrix<int_8> rtmx : " << endl;
    rtmx.Print(cout, -1); 
  }
  }
  td8.MinMax(min8, max8);
  if (PrtLev>2) {
    cout << " SymmetricMatrix<int_8> td8=tmx-tmx2 : " << endl;
    td8.Print(cout, -1); 
  }
  if ((min8!=0)||(max8!=0)) {
    cout << " !!! Test_2 ERROR : diff::min=" << min8 << " diff::max=" << max8 << endl;
    rc += 4;
  }
  else cout << " ---> Test_2 OK " << endl;


//----- Test 3
  cout << "3/ Test 3 : ConvertToStdmatrix() + operations on SymmetricMatrix<r_8>" << endl;
  SymmetricMatrix<r_8>  tdmx(NR);
  tdmx = RandomSequence(RandomSequence::Gaussian);
  TMatrix<r_8> dmx = tdmx.ConvertToStdMatrix();
  int nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (tdmx(r,c) != dmx(r,c))  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_3.a ERROR : nerr = " << nerr << endl;
    rc += 8;
  }
  else cout << " ---> Test_3.a OK " << endl;


  tdmx *= M_PI;    tdmx -= 0.67;
  dmx *= M_PI;    dmx -= 0.67;
  SymmetricMatrix<r_8>  tdmx2(dmx);
  SymmetricMatrix<r_8> tdd8 = tdmx-tdmx2;
  r_8 mind, maxd;
  tdd8.MinMax(mind, maxd);
  if (PrtLev>2) {
    cout << " SymmetricMatrix<r_8> tdmx : " << endl;
    tdmx.Print(cout, -1); 
    cout << " SymmetricMatrix<int_8> tdd8=tdmx-tdmx2 : " << endl;
    tdd8.Print(cout, -1); 
  }
  if ((fabs(mind)>1.e-19)||(fabs(maxd)>1e-19)) {
    cout << " !!! Test_3.b ERROR : diff::min=" << mind << " diff::max=" << maxd << endl;
    rc += 16;
  }
  else cout << " ---> Test_3.b OK " << endl;

  SymmetricMatrix<r_8>  tdmxa(NR);
  tdmxa = RegularSequence(2.,3.);
  TMatrix<r_8> dmxa = tdmxa.ConvertToStdMatrix();
  SymmetricMatrix<r_8>  tdmxb = (74.32*((5.89*tdmx-tdmxa)&&tdmxa)+14.5)/tdmx;
  TMatrix<r_8>  dmxb = (74.32*((5.89*dmx-dmxa)&&dmxa)+14.5)/dmx;
  SymmetricMatrix<r_8>  tdmxc(dmxb);
  tdd8 = tdmxb-tdmxc;
  tdd8.MinMax(mind, maxd);
  if (PrtLev>2) {
    cout << " SymmetricMatrix<r_8> tdmxb : " << endl;
    tdmx.Print(cout, -1); 
    cout << " SymmetricMatrix<int_8> tdd8=tdmxb-tdmxc : " << endl;
    tdd8.Print(cout, -1); 
  }
  if ((fabs(mind)>1.e-19)||(fabs(maxd)>1e-19)) {
    cout << " !!! Test_3.c ERROR : diff::min=" << mind << " diff::max=" << maxd << endl;
    rc += 32;
  }
  else cout << " ---> Test_3.c OK " << endl;

//----- Test 4
  cout << "4/ Test 4 : Matrix multiplication on SymmetricMatrix<r_8>" << endl;
  SymmetricMatrix<r_8>  tmxa(NR), tmxb(NR);
  tmxa = RandomSequence(RandomSequence::Gaussian);
  tmxb = RandomSequence(RandomSequence::Gaussian,3.);
  TMatrix<r_8> mxd = tmxa * tmxb;
  TMatrix<r_8> mxa = tmxa.ConvertToStdMatrix();
  TMatrix<r_8> mxb = tmxb.ConvertToStdMatrix();
  TMatrix<r_8> mxc = mxa * mxb;

  nerr = 0;
  for (sa_size_t c=0; c<mxd.NCols(); c++) {
    for (sa_size_t r=0; r<mxd.NRows(); r++) {
      if (fabs(mxd(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.a ERROR : nerr = " << nerr << endl;
    rc += 64;
  }
  else cout << " ---> Test_4.a OK " << endl;

  mxb = RandomSequence(RandomSequence::Gaussian, 9.);
  mxd = tmxa*mxb;
  mxc = mxa*mxb;
  if (PrtLev>2) {
    cout << " SymmetricMatrix<r_8> tmxa : " << endl;
    tmxa.Print(cout, -1); 
    cout << " TMatrix<r_8> mxb : " << endl;
    mxb.Print(cout, -1); 
    cout << " TMatrix<r_8> mxd=tmxa*mxb : " << endl;
    mxd.Print(cout, -1); 
  }

  nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (fabs(mxd(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.b ERROR : nerr = " << nerr << endl;
    rc += 128;
  }
  else cout << " ---> Test_4.b OK " << endl;

  mxd = mxb*tmxa;
  mxc = mxb*mxa;
  nerr = 0;
  for (sa_size_t c=0; c<tmx.NCols(); c++) {
    for (sa_size_t r=0; r<tmx.NRows(); r++) {
      if (fabs(mxd(r,c)-mxc(r,c))>1e-15)  nerr++;
    }
  }
  if (nerr>0) {
    cout << " !!! Test_4.c ERROR : nerr = " << nerr << endl;
    rc += 256;
  }
  else cout << " ---> Test_4.c OK " << endl;


  return rc;
}