#include <stdlib.h>
#include <stdio.h>

#include <exception>
#include <string>

#ifdef USEVECSTL
#include <vector>
#endif

using namespace std;

// class MyException : public exception
class MyException 
{
public:
  MyException(const char * msg) { _msg = msg; }
  // virtual ~MyException() { }
  string Msg() { return(_msg); }
private:
string _msg;
};

template<class T> class Matrix
{
  T* data;
  int siz_x, siz_y, offset, step_x, step_y;
  int size;
  
public:
  Matrix (int sx, int sy, int step=0, int offset=0, bool fg=false);
  ~Matrix ();

  inline T operator [] (int k) const { return data[k]; }
  inline T& operator [] (int k) { return data[k]; }

  inline T operator () (int ix, int iy) const { return data[ix+iy*siz_x]; }
  inline T& operator () (int ix, int iy) { return data[ix+iy*siz_x]; }

  inline T  elem (int ix, int iy) const { return data[offset+ix*step_x+iy*step_y]; }
  inline T& elem (int ix, int iy) { return data[offset+ix*step_x+iy*step_y]; }

  inline T  elemCk (int ix, int iy) const
    { if ((ix < 0) || (ix >= siz_x) || (iy < 0) || (iy >= siz_y)) 
                 throw MyException("Matrix<T> Out of bound");
      return data[offset+ix*step_x+iy*step_y]; }

  inline T& elemCk (int ix, int iy) 
    { if ((ix < 0) || (ix >= siz_x) || (iy < 0) || (iy >= siz_y)) 
                 throw MyException("Matrix<T> Out of bound");
      return data[offset+ix*step_x+iy*step_y]; }

  inline int getSize () const { return size; }
  inline int getSizeX () const { return siz_x; }
  inline int getSizeY () const { return siz_y; }

  // Addition et multiplication en utilisant elem()
  Matrix<T> * Add(Matrix<T> & v1, Matrix<T> & v2);
  Matrix<T> * Mult(Matrix<T> & v1, Matrix<T> & v2);
  // Addition et multiplication en utilisant elemCk()
  Matrix<T> * AddCk(Matrix<T> & v1, Matrix<T> & v2);
  Matrix<T> * MultCk(Matrix<T> & v1, Matrix<T> & v2);
  // Addition et multiplication en utilisant l'operateur []
  Matrix<T> * AddO1(Matrix<T> & v1, Matrix<T> & v2);
  Matrix<T> * MultO1(Matrix<T> & v1, Matrix<T> & v2);
  // Addition et multiplication en utilisant l'operateur ()
  Matrix<T> * AddO2(Matrix<T> & v1, Matrix<T> & v2);
  Matrix<T> * MultO2(Matrix<T> & v1, Matrix<T> & v2);
};

template<class T> Matrix<T>::Matrix (int sx, int sy, int step, int off, bool fg)
{
  int k;
  int s = offset+sx*sy*step;
  if (s < 1) s = 1;
  size = s;
  siz_x = sx; 
  siz_y = sy;
  offset = off;
  step_x = step;
  step_y = step*siz_x;
  data = new T[size];
  if (!fg)  return;

  T * p = data;
  for(k=0; k<s; k++) p[k] = (T)0;
}

template<class T> Matrix<T>::~Matrix()
{
delete[] data;
}



template<class T> Matrix<T> * Matrix<T>::Add(Matrix<T> &v1, Matrix<T> &v2)
{
int i,j;

for (i = 0; i < siz_x; i++)
  for (j = 0; j < siz_y; j++)
  elem(i,j) = v1.elem(i,j)+v2.elem(i,j);
return (this);
}

template<class T> Matrix<T> * Matrix<T>::Mult(Matrix<T> &v1, Matrix<T> &v2)
{
int i,j;

for (i = 0; i < siz_x; i++)
  for (j = 0; j < siz_y; j++)
  elem(i,j) = v1.elem(i,j)*v2.elem(i,j);
return (this);
}


template<class T> Matrix<T> * Matrix<T>::AddCk(Matrix<T> &v1, Matrix<T> &v2)
{
int i,j;

for (i = 0; i < siz_x; i++)
  for (j = 0; j < siz_y; j++)
  elem(i,j) = v1.elemCk(i,j)+v2.elemCk(i,j);
return (this);
}

template<class T> Matrix<T> * Matrix<T>::MultCk(Matrix<T> &v1, Matrix<T> &v2)
{
int i,j;

for (i = 0; i < siz_x; i++)
  for (j = 0; j < siz_y; j++)
  elem(i,j) = v1.elemCk(i,j)*v2.elemCk(i,j);
return (this);
}

template<class T> Matrix<T> * Matrix<T>::AddO1(Matrix<T> &v1, Matrix<T> &v2)
{
int i;
 for (i = 0; i < size; i++) (*this)[i] = v1[i] + v2[i];
return (this);
}

template<class T> Matrix<T> * Matrix<T>::MultO1(Matrix<T> &v1, Matrix<T> &v2)
{
int i;
 for (i = 0; i < size; i++) (*this)[i] = v1[i] * v2[i];
return (this);
}

template<class T> Matrix<T> * Matrix<T>::AddO2(Matrix<T> &v1, Matrix<T> &v2)
{
int i,j;

for (i = 0; i < siz_x; i++)
  for (j = 0; j < siz_y; j++)  (*this)(i,j) = v1(i,j) + v2(i,j);
return (this);
}

template<class T> Matrix<T> * Matrix<T>::MultO2(Matrix<T> &v1, Matrix<T> &v2)
{
int i,j;

for (i = 0; i < siz_x; i++)
  for (j = 0; j < siz_y; j++)  (*this)(i,j) = v1(i,j) * v2(i,j);
return (this);
}


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

/* --------------------------------------------------------------------- */
/* --------------------------- Main Program ---------------------------- */
/* --------------------------------------------------------------------- */

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

int pos, N, M, Mx, My, Off, Step, OPT, OPE, i;


if (narg < 2) { 
  printf("\n Usage: matrix Type(=1,2,3 Int,Float,Double) Ope(=1...10) [N [Mx,My,...] ] \n");
  printf("Ope: 1=Create/Delete 2=1+FillVect \n");
  printf("Ope: 3=AddO1  4=MultO1 (Using operator [k]) \n");
  printf("Ope: 5=AddO2  6=MultO2 (Using operator (i,j)) \n");
  printf("Ope: 7=Add    8=Mult   (Using elem()) \n");
  printf("Ope: 9=AddCk  10=MultCk (Using elemCk() - with bound checking) \n");
  printf("N: Number of operations (def= 100) \n");
  printf("Mx,My,Step,Offset: Matrix size (def= 300,200,1,0) \n\n");
  exit(0);
}

 InitTim();

OPT = 1; OPE = 1;
if (narg > 1)  OPT = atoi(arg[1]);
if ( (OPT < 1) || (OPT > 3) )  OPT = 1;

if (narg > 2)  OPE = atoi(arg[2]);
if ( (OPE < 1) || (OPE > 10) )  OPE = 1;

N = 100;
if (narg > 3)  N = atoi(arg[3]);
if (N < 1) N = 1;
if (N > 100000) N = 100000;
Mx = 300; 
My = 200;
Off = 0;
Step = 1;
if (narg > 4)  sscanf(arg[4], "%d,%d,%d,%d", &Mx, &My, &Step, &Off);
if (Mx < 100) Mx = 100;
if (My < 100) My = 100;
if (Mx > 10000) Mx = 10000;
if (My > 10000) My = 10000;
if (Step < 1) Step = 1;
if (Step > 5) Step = 5;
if (Off < 0) Off = 0;
if (Off > 1000) Off = 1000;

M = Mx*My*Step;


printf(" MatrixC++ TestSpeed Typ=%d Ope=%d N=%d MSz=%d\n", OPT, OPE, N,M);
printf("    Matrix Size X= %d  Y = %d  Step= %d Offset= %d \n", Mx, My, Step, Off);
if (OPE < 3) { /* Test creation / destruction */
int fg = OPE-1;
bool fgf = (fg != 0) ? true : false;
printf("\n\n Test new/delete Matrix<T> - fg= %d ( <> 0 ---> FillVec) \n", fg); 
 switch (OPT)
   {
   case 1 :
     {
       Matrix<int> * v;
       printf("Test %d new/delete Matrix<int>[%d] \n",N,M);
       for(i=0; i<N; i++) {
	 v = new Matrix<int>(Mx, My, Step, Off, fgf);
	 delete v;
       }
     }
     break;
   case 2 :
     {
       Matrix<float> * v;
       printf("Test %d  new/delete Matrix<float>[%d] \n",N,M);
       for(i=0; i<N; i++) {
	 v = new Matrix<float>(Mx, My, Step, Off, fgf);
	 delete v;
       }
     }
     break;
   case 3 :
     {
       Matrix<double> * v;
       printf("Test %d  new/delete Matrix<double>[%d] \n",N,M);
       for(i=0; i<N; i++) {
	 v = new Matrix<double>(Mx, My, Step, Off, fgf);
	 delete v;
       }
     }
     break;
  }
  PrtTim("Fin New/Delete ");
  printf(" .......... Fin de MatrixC++ ........... \n");
  return (0);
}

// ---------- Test Addition, Multiplication -------------

switch (OPT)
  {
  case 1 :
    {
    printf("\n\n Test %d operations Matrix<int>[%d] \n",N,M);

    Matrix<int> *v1,*v2,*v3;
    v1 = new Matrix<int>(Mx, My, Step, Off);
    v2 = new Matrix<int>(Mx, My, Step, Off);
    v3 = new Matrix<int>(Mx, My, Step, Off);
    PrtTim("Fin_Creation ");

    for(pos=0; pos<M; pos++)  
      { (*v1)[pos] = random()%1000;
      (*v2)[pos] = random()%5000; }

    PrtTim("Fin remplissage ");

    if (OPE == 3) {
      for(pos=0; pos<N; pos++)  
        v3->AddO1(*v1, *v2);  
      PrtTim("Fin Addition AddO1 operator [k]");
    }
    else if (OPE == 4) {
      for(pos=0; pos<N; pos++)  
        v3->MultO1(*v1, *v2);  
      PrtTim("Fin Multiplication MultO1  operator [k]");
    }
    else if (OPE == 5) {
      for(pos=0; pos<N; pos++)  
        v3->AddO2(*v1, *v2);  
      PrtTim("Fin Addition AddO2 operator (i,j)");
    }
    else if (OPE == 6) {
      for(pos=0; pos<N; pos++)  
        v3->MultO2(*v1, *v2);  
      PrtTim("Fin Multiplication MultO2 operator (i,j)");
    }
    else if (OPE == 7) {
      for(pos=0; pos<N; pos++)  
        v3->Add(*v1, *v2);  
      PrtTim("Fin Addition Add");
    }
    else if (OPE == 8) {
      for(pos=0; pos<N; pos++)  
        v3->Mult(*v1, *v2);  
      PrtTim("Fin Multiplication Mult");
    }
    else if (OPE == 9) {
      for(pos=0; pos<N; pos++)  
        v3->AddCk(*v1, *v2);  
      PrtTim("Fin Addition AddCk");
    }
    else if (OPE == 10) {
      for(pos=0; pos<N; pos++)  
        v3->MultCk(*v1, *v2);  
      PrtTim("Fin Multiplication MultCk");
    }

    printf("Result[1.2] I1= %d %d  I2= %d %d  I3= %d %d \n",
          v1->elem(1,0),v1->elem(2,0), v2->elem(1,0),v2->elem(2,0),
          v3->elem(1,0),v3->elem(2,0));
    printf("ResAdd[991-2] I1= %d %d  I2= %d %d  I3= %d %d \n",
          v1->elem(991,0),v1->elem(992,0), v2->elem(991,0),v2->elem(992,0),
          v3->elem(991,0),v3->elem(992,0));
    }
    break;

  case 2 :
    {
    printf("\n\n Test %d operations Matrix<float>[%d] \n",N,M);
    Matrix<float> *v1,*v2,*v3;
    v1 = new Matrix<float>(Mx, My, Step, Off);
    v2 = new Matrix<float>(Mx, My, Step, Off);
    v3 = new Matrix<float>(Mx, My, Step, Off);
    PrtTim("Fin_Creation ");

    for(pos=0; pos<M; pos++)  
      { (*v1)[pos] = (float)(random()%1000)/250.;
        (*v2)[pos] =  (float)(random()%5000)/250.; }

    PrtTim("Fin remplissage ");

    if (OPE == 3) {
      for(pos=0; pos<N; pos++)  
        v3->AddO1(*v1, *v2);  
      PrtTim("Fin Addition AddO1 operator [k]");
    }
    else if (OPE == 4) {
      for(pos=0; pos<N; pos++)  
        v3->MultO1(*v1, *v2);  
      PrtTim("Fin Multiplication MultO1  operator [k]");
    }
    else if (OPE == 5) {
      for(pos=0; pos<N; pos++)  
        v3->AddO2(*v1, *v2);  
      PrtTim("Fin Addition AddO2 operator (i,j)");
    }
    else if (OPE == 6) {
      for(pos=0; pos<N; pos++)  
        v3->MultO2(*v1, *v2);  
      PrtTim("Fin Multiplication MultO2 operator (i,j)");
    }
    else if (OPE == 7) {
      for(pos=0; pos<N; pos++)  
        v3->Add(*v1, *v2);  
      PrtTim("Fin Addition Add");
    }
    else if (OPE == 8) {
      for(pos=0; pos<N; pos++)  
        v3->Mult(*v1, *v2);  
      PrtTim("Fin Multiplication Mult");
    }
    else if (OPE == 9) {
      for(pos=0; pos<N; pos++)  
        v3->AddCk(*v1, *v2);  
      PrtTim("Fin Addition AddCk");
    }
    else if (OPE == 10) {
      for(pos=0; pos<N; pos++)  
        v3->MultCk(*v1, *v2);  
      PrtTim("Fin Multiplication MultCk");
    }

    printf("Result[1.2] F1= %g %g  F2= %g %g  F3= %g %g \n",
          v1->elem(1,0),v1->elem(2,0), v2->elem(1,0),v2->elem(2,0),
          v3->elem(1,0),v3->elem(2,0));
    printf("Result[991-2] F1= %g %g  F2= %g %g  F3= %g %g \n",
          v1->elem(991,0),v1->elem(992,0), v2->elem(991,0),v2->elem(992,0),
          v3->elem(991,0),v3->elem(992,0));
    }
    break;

  case 3 :
    {
    printf("\n\n Test %d operations Matrix<double>[%d] \n",N,M);
    Matrix<double> *v1,*v2,*v3;
    v1 = new Matrix<double>(Mx, My, Step, Off);
    v2 = new Matrix<double>(Mx, My, Step, Off);
    v3 = new Matrix<double>(Mx, My, Step, Off);
    PrtTim("Fin_Creation ");


    if (OPE == 3) {
      for(pos=0; pos<N; pos++)  
        v3->AddO1(*v1, *v2);  
      PrtTim("Fin Addition AddO1 operator [k]");
    }
    else if (OPE == 4) {
      for(pos=0; pos<N; pos++)  
        v3->MultO1(*v1, *v2);  
      PrtTim("Fin Multiplication MultO1  operator [k]");
    }
    else if (OPE == 5) {
      for(pos=0; pos<N; pos++)  
        v3->AddO2(*v1, *v2);  
      PrtTim("Fin Addition AddO2 operator (i,j)");
    }
    else if (OPE == 6) {
      for(pos=0; pos<N; pos++)  
        v3->MultO2(*v1, *v2);  
      PrtTim("Fin Multiplication MultO2 operator (i,j)");
    }
    else if (OPE == 7) {
      for(pos=0; pos<N; pos++)  
        v3->Add(*v1, *v2);  
      PrtTim("Fin Addition Add");
    }
    else if (OPE == 8) {
      for(pos=0; pos<N; pos++)  
        v3->Mult(*v1, *v2);  
      PrtTim("Fin Multiplication Mult");
    }
    else if (OPE == 9) {
      for(pos=0; pos<N; pos++)  
        v3->AddCk(*v1, *v2);  
      PrtTim("Fin Addition AddCk");
    }
    else if (OPE == 10) {
      for(pos=0; pos<N; pos++)  
        v3->MultCk(*v1, *v2);  
      PrtTim("Fin Multiplication MultCk");
    }


    printf("Result[1.2] D1= %g %g  D2= %g %g  D3= %g %g \n",
          v1->elem(1,0),v1->elem(2,0), v2->elem(1,0),v2->elem(2,0),
          v3->elem(1,0),v3->elem(2,0));
    printf("Result[991-2] D1= %g %g  D2= %g %g  D3= %g %g \n",
          v1->elem(991,0),v1->elem(992,0), v2->elem(991,0),v2->elem(992,0),
          v3->elem(991,0),v3->elem(992,0));
    }
    break;

  default:
    puts("Erreur d'option !");
    break;
  }


  
PrtTim("Fin de Matrix/C++ ");
printf(" .......... Fin de MatrixC++ ........... \n");
return (0);
}

/*
#ifdef DECCXX
#pragma define_template Matrix<int>
#pragma define_template Matrix<float>
#pragma define_template Matrix<double>
#endif

#if defined(GNUGCC) || defined (HPaCC)
template class Matrix<int>;
template class Matrix<float>;
template class Matrix<double>;
#endif

*/