source: Sophya/trunk/SophyaProg/Tests/arrt.cc@ 1555

#include "machdefs.h"

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

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

static void tstarr(bool tio);
static void tstioarr(TArray<int_4> & ia, TArray<int_4> & ib, TArray<int_4> & ic);
static void tstmtx(int n);
void _ZZ_TestTMatrixRC_YY_(TMatrix<r_8> & m);

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

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


  if (narg < 2) {
    cout << " arrt - TArray<T> Test - Usage arrt a/w/m/z [size=5] [prtlev=1] [nprt=50]\n" 
         << " a: Simple Array Test  w: a+ PPersist Test (FIO_TArray<T>) \n" 
         << " m: Matrix and Vector Test \n" 
         << " z: Appel_ZZ_TestTMatrixRC_YY_ " << endl;
    exit(0);
  }

  int opt = *(arg[1]);
  int n = 5;
  if (narg > 2) n = atoi(arg[2]);
  uint_4 prtlev = 1;
  uint_4 nprt = 50;
  if (narg > 3) prtlev = atoi(arg[3]);
  if (narg > 4) nprt = atoi(arg[4]);

  BaseArray::SetMaxPrint(nprt, prtlev);

  try {
    if ((opt == 'w') || (opt == 'a')) tstarr((opt=='w'));
    else if (opt == 'z') {
      if (n<5) n = 5;
      TMatrix<r_8> tza(n,n+3);
      tza = RegularSequence(2.,1.5);
      _ZZ_TestTMatrixRC_YY_(tza);
      TMatrix<r_8> tza2 = tza(Range(1,3), Range(2,5));
      _ZZ_TestTMatrixRC_YY_(tza2);
    }
    else tstmtx(n);
  }
  catch (PThrowable exc) {
    cerr << " catched Exception " << exc.Msg() << endl;
  }  
  catch (...) {
    cerr << " catched unknown (...) exception " << endl; 
  }  
  PrtTim("--- End of arrt ---");
  cout << " ---------------  END of ArrayTest programme -------------- " << endl;
}


void tstmtx(int n)
{
  TMatrix<r_4> ma(4,7);
  ma = RegularSequence(10., 2.);
  cout << "\n ::::::  Testing TMatrix - TVector ::::::  " << endl;  
  TMatrix<r_4> sma1 = ma(Range(1,2), Range(2,4));
  cout << " ----- Matrix ma = \n " << ma << endl;
  sma1 = 3.14;
  cout << " ----- Matrix ma (Avec sma1=ma(Range(1,2), Range(2,3))=3.14) = \n " << ma << endl;
  TMatrix<r_4> mb(2,3);
  mb = 50.;
  TMatrix<r_4> mc = mb+sma1;
  cout << " ----- Matrix mc = mb(=50.)+sma1(=3.14) = " << mc << endl;
  TMatrix<r_4> md = (float)100+mc;
  cout << " ----- Matrix md=100.+mc = " << md << endl;
  cout << " ----- Matrix md-50 = " << md-(float)50. << endl;

  TMatrix<int_4> ima(3,10);
  ima = RegularSequence(0.01,1.);
  cout << " ----- Matrix ima(3,10) = RegularSequence(0.01,1.) = " << ima << endl;
  cout << " ----- Matrix  ima.Column(3) = \n" << ima.Column(3) << endl;
  TVector<int_4> vcol = ima.Column(7);
  cout << " ----- Vector  vcol=ima.Column(7) = \n" << vcol << endl;
  cout << " ----- Matrix  ima.Row(0) = \n" << ima.Row(0) << endl;
  TVector<int_4> vlig;
  vlig = ima.Row(1);
  cout << " ----- Vector  vlig=ima.Row(1) = \n" << vlig << endl;

  TMatrix<int_4> imb(3,10);
  imb = RegularSequence(100.01,10.);
  cout << " ----- Matrix imb(3,10) = RegularSequence(100.01,10.) = " << imb << endl;
  TMatrix<int_4> imc = ima+imb;
  cout << " ----- Matrix imc=ima+imb = " << imc << endl;
  TMatrix<int_4> imd(3,10);
  TMatrix<int_4> imd2(3,10);
  imd += 25;
  imd *= 2;
  imd2 = 2;
  imd.MulElt(imd2);
  cout << " ----- Matrix imd (? = 100) " << imd << endl;
  imb -= imd;
  imc -= imd;
  cout << " ----- Matrix ((imc-=imd)-(imb-=imd))-ima (== 0 ?) " 
       << (imc-imb)-ima << endl;
  TVector<r_4> va(16, BaseArray::RowVector);
  va = RegularSequence(0.,1.);
  cout << " ----- Vector va = RegularSequence(0.,1.) = " << va << endl;
  va(Range(2,0,6,2)) *= 10.;
  cout << " ----- Vector va apres va(Range(2,0,6,2)) *= 10. " << va << endl;
  va(Range(0,5)) += 100.;
  cout << " ----- Vector va apres va(Range(0,6)) += 100. " << va << endl;

  TMatrix<r_4> esm(3,5);
  EnumeratedSequence ess;
              ess ,
         11 , 12 , 13, 14, 15,  
         21 , 22 , 23, 24, 25,
         31 , 32 , 33, 34, 35  ;  
  esm = ess;
  cout << "\n TMatrix<r_4> esm(3,5) = EnumeratedSequence(), 11 , 12 , 13, 14, 15 ... " 
       << esm << endl;
  r_4 emin,emax,esum,eprod,esumx2;
  esm.MinMax(emin, emax);  
  esum = esm.Sum();
  eprod = esm.Product();
  esumx2 = esm.SumX2();

  cout << " Min= " << emin << "  Max= " << emax 
       << " Sum= " << esum << " Prod=" << eprod << " SumX2=" << esumx2 << endl;

  cout << "\n :::::: Testing Matrix operations :::::: " << endl;
  TMatrix<r_4> ra(2,2), ira(2,2);
  /*  ra(0,0) = 2.;   ra(0,1) = 1.;
  ra(1,0) = 4.;   ra(1,1) = 3.;
  */
  ira(0,0) = 1.5;   ira(0,1) = -0.5;
  ira(1,0) = -2.;    ira(1,1) = 1.;
  EnumeratedSequence es;
  es, 2 , 1 , 4 , 3 ;
  ra = es ;
  
  cout << " ----- Matrix ra = \n " << ra << endl; 
  cout << " ----- Matrix ira = \n " << ira << endl; 
  cout << " ----- Matrix ira*ra = (= Identity ?) " << ira*ra << endl; 


  cout << "\n ::::::: Testing Matrix Inversion ::::::: " << endl;
  SimpleMatrixOperation<r_4> smr;
  TMatrix<r_4> cira = smr.Inverse(ra);
  cout << " ----- Matrix Inverse(ra) =  " << cira << endl;
  cout << " ----- ira-Inverse(ra) = " << ira-cira << endl;

  Matrix a(n,n);

    
  //  for(int i=0; i<n; i++)
  //    for(int j=0; j<n; j++)  a(j,i) = GauRnd(0., 1.);
  a = RandomSequence(RandomSequence::Gaussian, 0., 2.5);
  cout << " ----- Matrix a=RandomSequence(RandomSequence::Gaussian, 0., 2.5) = \n " << a << endl;
  Vector x(n);
  x = RegularSequence(1.,3.);
  Vector b = a*x;
  cout << " ----- Vector x = \n " << x << endl;
  cout << " ----- Vector b = a*x = \n " << b << endl;
  SimpleMatrixOperation<r_8> smo;
  Matrix inva = smo.Inverse(a);
  cout << " ----- Matrix Inverse(a) = \n " << inva << endl;
  cout << " ----- Matrix a*Inverse(a) = \n " << inva*a << endl;
  cout << " ----- Matrix Inverse(a)*b (=Inv(a)*a*x) = \n " << inva*b << endl;
  cout << " ----- Matrix x-Inverse(a)*b = (=0 ?)\n " << x-inva*b << endl;
}

void tstarr(bool tio)
{
  cout << "\n -----> Testing TArray <---- " << endl;  
  // We create a integer array SizeX=7, SizeY=5
  TArray<int_4> ia(7,5); 
  // We fill it with a sequence of numbers starting at 10., with step = 2.
  ia = RegularSequence(10., 2.);
  cout << " ----- Array IA = \n " << ia << endl;
  // sub array extraction, Range(2,4) : starting position=2 , End=4
  TArray<int_4> ic = ia(Range(2,3),Range(1,3),Range(0));
  cout << " ----- Array IC IA(Range(2,3),Range(1,3)) = \n " << ic << endl;
  // we set the sub-array to zero, this should reflect in the original array
  // sub-arrays share their data with parent array
  ic = 0;
  cout << " ----- Array IC Apres (=0) = \n " << ic << endl;
  cout << " ----- Array IA Apres IC=0 = \n " << ia << endl;

  cout << "  :::: 3 Dim arrays ::::: " << endl;
  TArray<int_4>::SetMaxPrint(1000);
  // Creating 3-dim array (X=8 x Y=7 x Z=2) , filling it with 5
  TArray<int_4> ib(8,7,2);
  ib = 5;
  cout << " ----- Array IB = \n " << ib << endl;
  // Sub array extraction  X from 1 , size 4 - Y from 2 , size 3 , in Z default, from 0, size 1
  // we multiply this sub-array elements by 3
  ib(Range(1,4),Range(2,4), Range(0)) *= 3;
  cout << " -- Array IB , Apres ib(Range(1,4),Range(2,4))*=3 : " << endl;
  cout << ib;
  cout << " -- Array ib(Range(1,4),Range(2,4)) =  : " 
       << ib(Range(1,4),Range(2,4), Range(0)) << endl;

  // Creating a double array X=5 x Y=2
  TArray<r_4> fa(5,2);
  // fill it up with a sequence of 0. to 1.
  fa = RegularSequence(0.,1./(5*2));
  cout << " ------ TArray<r_4> fa(5,2) = \n" << fa << endl;
  // Create a new array from the original array , multiplying it by 2*Pi
  TArray<r_4> fa2 = fa*(float)(2.*M_PI); 
  cout << " ------ TArray<r_4> fa2=fa*2*Pi = \n" << fa2 << endl;
  // Compute sin(fa2) cos(fa2)
  cout << " ------ sin(fa2=fa*2*Pi) = \n" << Sin(fa2) << endl;
  cout << " ------ cos(fa2=fa*2*Pi) = \n" << Cos(fa2) << endl;

  //  Test des fonctions sur complexes 
  TVector<r_4> p_real(10, BaseArray::RowVector), p_imag(10, BaseArray::RowVector);
  p_real = RegularSequence(0., 0.5);
  p_imag = RegularSequence(0., 0.25);
  TVector< complex<r_4> > zvec = ComplexArray(p_real, p_imag);
  cout << " ------ zvec = ComplexArray(p_real, p_imag) -------- "  << endl;
  cout << " :: p_real= " << p_real; 
  cout << " :: p_imag= " << p_imag;
  cout << " :: zvec= " << zvec;
  cout << " :: real(zvec) = " << real(zvec) ;
  cout << " :::: imag(zvec) = " << imag(zvec) ;
  cout << " :::: module2(zvec) = " << module2(zvec) ;
  cout << " :::: module(zvec) = " << module(zvec) ;
  cout << " :::: phase(zvec) = " << phase(zvec) ;
  cout << " ----------------------------------------------------- " << endl;

  if (tio) tstioarr(ia, ib, ic);
}

void tstioarr(TArray<int_4> & ia, TArray<int_4> & ib, TArray<int_4> & ic)
{
  cout << " ------ tstioarr(TArray<int_4> & ia, TArray<int_4> & ib, TArray<int_4> & ic) ---- " << endl;
  {
  cout << " >>>>>> Writing in  arrt.ppf <<<<<<< " << endl;
  POutPersist pos("arrt.ppf");
  // We write  the three arrays in the stream 
  pos << ia << ib << ic;
  cout << " >>>>>> Writing in  arrtn.ppf with names <<<<<<<" << endl;
  POutPersist posn("arrtn.ppf");
  string tag = "ArrIA";
  posn.PutObject(ia, tag);
  tag = "ArrIB";
  posn.PutObject(ib, tag);
  tag = "ArrIC";
  posn.PutObject(ic, tag);
  }

  {
  cout << " >>>>>>> Reading from  arrt.ppf  <<<<< " << endl;
  PInPersist pis("arrt.ppf");
  TArray<int_4> iaa, ibb, icc; 
  // We read the three arrays from the stream 
  pis >> iaa >> ibb >> icc;
  cout << " ----- Array IAA = \n " << iaa << endl;
  cout << " ----- Array IBB = \n " << ibb << endl;
  cout << " ----- Array ICC = \n " << icc << endl;
  icc = 12;
  cout << " ----- Array ICC (=12) = \n " << icc << endl;
  cout << " ----- Array IAA (ICC=12) = \n " << iaa << endl;
  }

  {
  cout << " >>>>>>> Reading from  arrtn.ppf  <<<<< " << endl;
  PInPersist pis("arrtn.ppf");
  TArray<int_4> iaa, ibb, icc; 
  // We read the three arrays from the stream 
  string tag = "ArrIC";
  pis.GetObject(icc, tag);
  tag = "ArrIB";
  pis.GetObject(ibb, tag);
  tag = "ArrIA";
  pis.GetObject(iaa, tag);

  cout << " ----- Array IAAA = \n " << iaa << endl;
  cout << " ----- Array IBBB = \n " << ibb << endl;
  cout << " ----- Array ICCC = \n " << icc << endl;
  icc = 68;
  cout << " ----- Array ICCC (=12) = \n " << icc << endl;
  cout << " ----- Array IAAA (ICC=12) = \n " << iaa << endl;
  }


}