#include<iostream>
#include"CLHEP/Matrix/Matrix.h"
#include"CLHEP/Matrix/SymMatrix.h"

/* -------------------------------------------------
    Test comparatif (pour SOPHYA) des matrices CLHep
                            R. Ansari , Juin 2005

    Pour compiler :
    definir flag/chemin include et librariries CLHep,
    Au CC-IN2P3 par exemple, 
    csh> set inchep = '-I/usr/local/lhcxx/CLHEP/pro/include/'
    csh> set libhep = '-L/usr/local/lhcxx/CLHEP/pro/lib -lCLHEP'
    csh> cc -c timing.c 
    csh> g++ -c $inchep mtxclhep.cc
    csh> g++ -o mtxclhep mtxclhep.o timing.o $libhep

    Execution :
    Acces aux elements et operation add/multcst
    csh> time ./mtxclhep 50 1000 500
    Inversion / multiplication
    csh> time ./mtxclhep 1 1000 1000
    ------------------------------------------- */


extern "C" {
void InitTim(void);
void PrtTim(const char *Comm);
}

using namespace std;

int main(int narg, char* arg[])
{
  
  if (narg < 3) {
    cout << "\n mtxclhep - missing arguments\n" 
	 << "     Usage:  mtxclhep NLoop NRow NCol \n"  
	 << " if nrow==ncol ==>  test inverse/mtx multiply (one time)" << 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 << " ---- mtxclhep: NLoop=" << N << " NRow= " << nrow 
       << " NCol= " << ncol << " ---- " << endl;
  
  InitTim();
  srand48(10);

  //ATTENTION :  L'indexage commence a 1 
  try {
    int i,j,k;

    HepMatrix m1(nrow, ncol);
    HepMatrix m2(nrow, ncol);
    HepMatrix m3(nrow, ncol);
    HepMatrix m4(nrow, ncol);
    HepMatrix m5(nrow, ncol);

 
    for(k=0; k<N; k++) {
      double xxr = drand48();
      double yyr = 2.*drand48();
      for(i=1; i<=nrow; i++)
	for(j=1; 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) -- ");
 
    for(k=0; k<N; k++) {
      double c1 = drand48() + 1.2;
      double c2 = drand48() + 3.5;
      double c3 = drand48() + 6.7;
      m5 = m1*c1+m2*c2+m3*c3;
    }
    cout << " (3) Add/CstMult: m1*c1 + m2*c2 + m3*c3 done ---" << endl;
    PrtTim(" (3) -- ");

  
  
    if (nrow == ncol) {
      cout << " mtxclhep/inverse , mtx multiply A(N,N)  N=" << nrow << endl;
      HepMatrix A(nrow, nrow);
      
      for(i=1; i<=nrow; i++)
	for(j=0; j<=nrow; j++)   A(i,j) = drand48();
      
      cout << " (91)  A(i,j) = drand48() done ---" << endl;
      PrtTim(" (91) -- ");
      
      int ierr;
      HepMatrix Ainv = A.inverse(ierr);
      cout << " (92)  Ainv = A.inverse() done ---" << endl;
      PrtTim(" (92) -- ");
      if (ierr == 0) {
	cout << " A.inverse(ierr) OK computing A*Ainv ..." << endl; 
	HepMatrix mxx = A*Ainv;
	cout << " (93)  mxx = A*Ainv done ---" << endl;
	PrtTim(" (93) -- ");
      }
    }
  }
  catch (std::exception & sex) {
    cerr << " mtxclhep/ catched std::exception what()=" << sex.what() << endl;
  }  
  catch (...) {
    cerr << " mtxclhep/catched unknown (...) exception " << endl; 
  }  
  cout << " ------ FIN mtxclhep -------- " << endl;
  return 0;
}