#include "machdefs.h"

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

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

#include "pexceptions.h"
#include "srandgen.h"

#include "sophyainit.h"
#include "ctimer.h"
#include "timing.h"

#include "smtx.h"

//--------------------------------------------------------
//   Test d'expression template sur tableaux 
//       Reza - Avril 2003
//       revisite : Mars 2012 
//--------------------------------------------------------

//-------  Le code de smtx.cc

#include "smtx.cc"

//-------
// ------------------------------------------------
//   programme de test template expressions 
// ------------------------------------------------

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

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

  if (narg < 3) {
    cout << "\n te/ TemplateExpression Test - missing arguments\n" 
	 << "     Usage: te NLoop NRow NCol \n"  << endl;
    return 1;
  }
  int N = 10;
  int nrow = 2; 
  int ncol = 2;

  N = atoi(arg[1]);

  nrow = atoi(arg[2]);
  ncol = atoi(arg[3]);

  cout << " te/ TemplateExpression Test - NLoop = " << N 
       << "  NRow=" << nrow << " NCol= " << ncol 
       << "  Size=nr*nc= " << nrow*ncol << endl;

  try {
    cout << " ---- Verification operation ------ " << endl;
    SimpleMatrix<int_4> mx1(4,6);
    mx1 = 3;
    mx1 *= 2;
    cout << " mx1=3; mx1*=2; : " << endl;
    cout << mx1 << endl;
    SimpleMatrix<int_4> mx2(4,6);
    mx2.Init(1,1);
    SimpleMatrix<int_4> mx3(4,6), mx4;
    mx3.Init(1000,100);
    mx4 = mx3+mx2*2 ;
    cout << " mx2, mx3, mx4=mx3+2mx2 " << endl;
    cout << mx2 << endl << mx3 << endl << mx4 << endl;

    cout << ">> Test with SimpleMatrix<T>/JET " << endl;
    PrtTim(" (0) -- ");

    int i,j,k;
    {
    SimpleMatrix<r_8> m1(nrow, ncol);
    SimpleMatrix<r_8> m2(nrow, ncol);
    SimpleMatrix<r_8> m3(nrow, ncol);
    SimpleMatrix<r_8> m4(nrow, ncol);
    SimpleMatrix<r_8> m5(nrow, ncol);
    for(k=0; k<N; k++) {
      double xxr = drand01();
      double yyr = 2.*drand01();
      for(i=0; i<nrow; i++)
	for(j=0; j<ncol; j++) {
	  m1(i,j) = k*300+10.*i+j+xxr;	
	  m2(i,j) = k*550+20.*i+2.*j+yyr;
	  m3(i,j) = k*860.+40.*i+7.*j+yyr*3.14;
	}
    }

    cout << " (1) Element access done ---" << endl;
    PrtTim(" (1) -- ");
    // Calcul m1*c1 + m2*c2 + m3*c3;
    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
      m5 = m1;  m5.MulCst(c1); 
      m4 = m2;  m4.MulCst(c2);
      m5.AddElt(m4);
      m4 = m3;  m3.MulCst(c2);
      m5.AddElt(m4);
    }
    cout << " (2) Add/MulElt/Cst m1*c1 + m2*c2 + m3*c3 done ---" << endl;
    PrtTim(" (2) -- ");

    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
      m5 = m1*c1+m2*c2+m3*c3;
    }
    cout << " (3) JET: m1*c1 + m2*c2 + m3*c3 done ---" << endl;
    PrtTim(" (3) -- ");
    // Calcul m1*c1 + m1*(m2+c2) + m2*c3; 
    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
      m3 = m1;  m3.MulCst(c1); 
      m4 = m2;  m4.AddCst(c2); m4.MulElt(m1);
      m5 = m2;  m5.MulCst(c3);
      m5.AddElt(m3); m5.AddElt(m4);
    }
    cout << " (4) Add/MulElt/Cst m1*c1 + m1*(m2+c2) + m2*c3 done ---" << endl;
    PrtTim(" (4) --");

    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
      //      m4 = SMExprAdd< r_8, SMExprMtx<r_8>, SMExprMtx<r_8> >( SMExprMtx<r_8>(m1), SMExprMtx<r_8> (m2) );  
    // Calcul m1*c1 + m1*(m2+c2) + m2*c3; 
      m5 = m1*c1 + m1*(m2+c2) + m2*c3;
    }
    cout << " (5) JET: m5 = m1*c1 + m1*(m2+c2) + m2*c3 ---" << endl;
    PrtTim(" (5) -- ");

    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
    // Calcul m1*c1 + m1*(m2+c2) + m2*c3; 
      m5 = m1*c1 + m2 + c2*Sin(c3*m3+m1);
    }
    cout << " (6) JET: m5 = m1*c1 + m2 + c2*Sin(c3*m3+m1) ---" << endl;
    PrtTim(" (6) -- ");

    }

    
    // Test avec les TArray / TMatrix
    {
    cout << " >>>>> test with SOPHYA::TArray / TMatrix " << endl;
    TMatrix<r_8> m1(nrow, ncol);
    TMatrix<r_8> m2(nrow, ncol);
    TMatrix<r_8> m3(nrow, ncol);
    TMatrix<r_8> m4(nrow, ncol);
    TMatrix<r_8> m5(nrow, ncol);
    for(k=0; k<N; k++) {
      double xxr = drand01();
      double yyr = 2.*drand01();
      for(i=0; i<nrow; i++)
	for(j=0; j<ncol; j++) {
	  m1(i,j) = k*300+10.*i+j+xxr;	
	  m2(i,j) = k*550+20.*i+2.*j+yyr;
	  m3(i,j) = k*860.+40.*i+7.*j+yyr*3.14;
	}
    }
    cout << " (21) Element access done ---" << endl;
    PrtTim(" (21) -- ");
    // Calcul m1*c1 + m2*c2 + m3*c3;
    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
      /* Changement interface SOPHYA  : 2004-2005  
	 m5 = m1;  m5.Mul(c1); 
	 m4 = m2;  m4.Mul(c2);
	 m5.AddElt(m4);
	 m4 = m3;  m3.Mul(c3);
	 m5.AddElt(m4);
      */
      m1.MulCst(c1, m5); 
      m2.MulCst(c2, m4);
      m5.AddElt(m4, m5);
      m3.MulCst(c3,m4);
      m4.AddElt(m5,m4);
    }
    cout << " (22) Add/MulElt/Cst m1*c1 + m2*c2 + m3*c3 done ---" << endl;
    PrtTim(" (22) -- ");

    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
      m5 = m1*c1+m2*c2+m3*c3;
    }
    cout << " (23) m1*c1 + m2*c2 + m3*c3 done ---" << endl;
    PrtTim(" (23) -- ");

    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
    // Calcul m1*c1 + m1*(m2+c2) + m2*c3; 
      m5 = m1*c1 + m1 && (m2+c2) + m2*c3;
    }
    cout << " (25) m5 = m1*c1 + m1*(m2+c2) + m2*c3 ---" << endl;
    PrtTim(" (25) -- ");

    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
      m5 = m1*c1 + m2 + c2*Sin(c3*m3+m1);
    }
    cout << " (26) m5 = m1*c1 + m2 + c2*Sin(c3*m3+m1) ---" << endl;
    PrtTim(" (26) -- ");

    }

    // Test avec des double * x;
    {
    cout << " >>>>> test with r_8 * p = new double[size] " << endl;
    r_8 * m1 = new r_8[nrow*ncol];
    r_8 * m2 = new r_8[nrow*ncol];
    r_8 * m3 = new r_8[nrow*ncol];
    r_8 * m5 = new r_8[nrow*ncol];
    for(k=0; k<N; k++) {
      double xxr = drand01();
      double yyr = 2.*drand01();
      for(i=0; i<nrow*ncol; i++) {
	  m1[i] = k*300+10.*i+j+xxr;	
	  m2[i] = k*550+20.*i+2.*j+yyr;
	  m3[i] = k*860.+40.*i+7.*j+yyr*3.14;
      }
    }
    cout << " (31) p[i] Element access done ---" << endl;
    PrtTim(" (31) -- ");
    // Calcul m1*c1 + m2*c2 + m3*c3;
    for(k=0; k<N; k++) {
      double c1 = drand01() + 1.2;
      double c2 = drand01() + 3.5;
      double c3 = drand01() + 6.7;
      for(i=0; i<nrow*ncol; i++) 
	m5[i] = c1*m1[i]+c2*m2[i]+c3*m3[i];
    }
    cout << " (32) p[i] m1*c1 + m2*c2 + m3*c3 done ---" << endl;
    PrtTim(" (32) -- ");

    delete[] m1;
    delete[] m2;
    delete[] m3;
    delete[] m5;
    }
  }
  catch (PThrowable exc) {
    cerr << " te/catched Exception " << exc.Msg() << endl;
  }  
  catch (...) {
    cerr << " te/catched unknown (...) exception " << endl; 
  }  

  cout << "\n --------------------------------------------------------" << endl;
  PrtTim("--- End of te ---");
  cout << " ----  END of te/TemplateExpression programme ------ " << endl;
  return 0;
}