// matrix.hh

#include <iostream>

#ifndef		REALMATRIX_HH_
# define	REALMATRIX_HH_

// la class TYPE ne peut etre que de type nombre (ex: double, int, long,...)

template <class TYPE> class Matrix
{
public:
  // Constructor
  Matrix(const unsigned int i, const unsigned int j)
  {
    this->_i = i;
    this->_j = j;
    this->_lines = new _line[i];	// Create lines...
    for (unsigned int k = 0; k < i; ++k)
      this->_lines[k] = new TYPE[j];	// ...and create column
    for (unsigned int k = 0; k < i; ++k)
      for (unsigned int l = 0; l < j; ++l)
	this->_lines[k][l] = 0;
  }

  // Constructor by copy
  Matrix(const Matrix& mat)
  {
    this->_i = mat._i;
    this->_j = mat._j;
    this->_lines = new _line[this->_i];
    for (unsigned int k = 0; k < this->_i; ++k)
      {
	this->_lines[k] = new TYPE[this->_j];
	for (unsigned int l = 0; l < this->_j; ++l)
	  this->_lines[k][l] = mat._lines[k][l];
      }
  }

  // Destructor
  ~Matrix()
  {
    if (this->_i != 0 || this->_j != 0)
      {
	for (unsigned int i = 0; i < this->_i; i++)
	  delete [] this->_lines[i];
	delete [] this->_lines;
      }
  }

  // Affected operator
  Matrix& operator=(const Matrix& mat)
  {
    if (&mat != this)
      {
	if (mat._i != this->_i || mat._j != this->_j)	// Test dim
	  {
	    this->~Matrix();	// Destroy...
	    this->_i = mat._i;
	    this->_j = mat._j;
	    this->_lines = new _line[this->_i];	// ...and realloc
	    for (unsigned int i = 0; i < this->_i; ++i)
	      this->_lines[i] = new TYPE[this->_j];
	  }
	for (unsigned int i = 0; i < this->_i; ++i)	// Copy
	  for (unsigned int j = 0; j < this->_j; ++j)
	    this->_lines[i][j] = mat._lines[i][j];
      }
    return (*this);
  }
  
  // Add elements to matrix
  Matrix&	add_elems(const int& nb)
  {
    Matrix mat(this->_i + nb, this->_j + nb);
    for (unsigned int i = 0; i < this->_i; i++)
      for (unsigned int j = 0; j < this->_j; j++)
	mat(i, j) = this->_lines[i][j];
    *this = mat;
    return (*this);
  }
  
  // Add matrix to matrix
  Matrix&	add_mat(const Matrix& mat)
  {
    unsigned int i = this->_i;
    unsigned int j = this->_j;
    this->add_elems(mat._i);
    for (unsigned int k = i; k < this->_i; k++)
      for (unsigned int l = j; l < this->_j; l++)
	this->_lines[k][l] = mat(k - i, l - j);
    return (*this);
  }

  // Add of matrix
  Matrix		operator+(const Matrix& mat)
  {
    Matrix mat_tmp(this->_i, mat._j);
    if (this->_i != mat._i || this->_j != mat._j)
      {
	std::cerr << "Error of addiction of 2 matrix (dimensions)"\
		  << std::endl;
	return (mat_tmp);
      }
    for (unsigned int i = 0; i < this->_i; ++i)
      for (unsigned int j = 0; j < this->_j; ++j)
	mat_tmp._lines[i][j] = this->_lines[i][j] + mat._lines[i][j];
    return (mat_tmp);
  }

  Matrix		operator+(const Matrix& mat) const
  {
    Matrix mat_tmp(this->_i, mat._j);
    if (this->_i != mat._i || this->_j != mat._j)
      {
	std::cerr << "Error of addiction of 2 matrix (dimensions)"\
		  << std::endl;
	return (mat_tmp);
      }
    for (unsigned int i = 0; i < this->_i; ++i)
      for (unsigned int j = 0; j < this->_j; ++j)
	mat_tmp._lines[i][j] = this->_lines[i][j] + mat._lines[i][j];
    return (mat_tmp);
  }

  // Sub matrix
  Matrix		operator-(const Matrix& mat)
  {
    Matrix mat_tmp(this->_i, mat._j);
    if (this->_i != mat._i || this->_j != mat._j)
      {
	std::cerr << "Error of substraction of 2 matrix (dimensions)"\
		  << std::endl;
	return (mat_tmp);
      }
    for (unsigned int i = 0; i < this->_i; ++i)
      for (unsigned int j = 0; j < this->_j; ++j)
	mat_tmp._lines[i][j] = this->_lines[i][j] - mat._lines[i][j];
    return (mat_tmp);
  } 

  Matrix		operator-(const Matrix& mat) const
  {
    Matrix mat_tmp(this->_i, mat._j);
    if (this->_i != this->_j || this->_i != mat._i || mat._i != mat._j)
      {
	std::cerr << "Error of substraction of 2 matrix (dimensions)"\
		  << std::endl;
	return (mat_tmp);
      }
    for (unsigned int i = 0; i < this->_i; ++i)
      for (unsigned int j = 0; j < this->_j; ++j)
	mat_tmp._lines[i][j] = this->_lines[i][j] - mat._lines[i][j];
    return (mat_tmp);
  }
  
  // Multiplication of matrix
  Matrix		operator*(const Matrix& mat)
  {
    Matrix mat_tmp(this->_i, mat._j);
    if (this->_j != mat._i)	// Check dimension
      {
	std::cerr << "Error of produce of 2 matrix (dimensions)"\
		  << std::endl;
	return (mat_tmp);
      }
    for (unsigned int i = 0; i < this->_i; ++i)
      {
	for (unsigned int j = 0; j < mat._j; ++j)
	  {
	    TYPE res = 0;
	    // Produce lines and column of matrix
	    for (unsigned int k = 0; k < this->_j; ++k)
	      res += this->_lines[i][k] * mat._lines[k][j];
	    mat_tmp._lines[i][j] = res;
	  }
      }
    return (mat_tmp);
  }

  // Multiplication with real
  Matrix&		operator*(const TYPE& nbr)
  {
    for (unsigned int i = 0; i < this->_i; i++)
      for (unsigned int j = 0; j < this->_j; j++)
	{
	  this->_lines[i][j] *= nbr;
	}
    return (*this);
  }

  // Divide with real
  Matrix&		operator/(const TYPE& nbr)
  {
    for (unsigned int i = 0; i < this->_i; i++)
      for (unsigned int j = 0; j < this->_j; j++)
	{
	  this->_lines[i][j] /= nbr;
	}
    return (*this);
  }  

  //Get dimension of matrix
  const int	get_dim() const
  {
    if (this->_i != this->_j)
      {
	std::cerr << "Matrix isn't nxn, cannot give dimension" << std::endl;
	return (0);
      }
    return (this->_i);
  }
  
  // Operator to change mat(i, j)
  // to can write mat(i, j) = ...
  TYPE&	operator()(unsigned int i, unsigned int j)
  {
    return (this->_lines[i][j]);
  }
  
  // Operator to access mat(i, j)
  // to can write int x = mat(i, j)
  TYPE	operator()(unsigned int i, unsigned int j) const
  {
    return (this->_lines[i][j]);
  }
  
  // Test matrix nul
  const bool	test_nul() const
  {
    for (unsigned int i = 0; i < this->_i; i++)
      for (unsigned int j = 0; j < this->_j; j++)
	if (this->_lines[i][j] != 0)
	  return (0);
    return (1);
  }
  
  // Test diagonal matrix nul
  const int	test_diag() const
  {
    for (unsigned int i = 0; i < this->_i; i++)
      if (this->_lines[i][i] != 0)
	return (i);
    return (-1);
  }
  
  // Test diagonal matrix nul with element
  const int	test_diag(const unsigned int& j) const
  {
    for (unsigned int i = 0; i < this->_i; i++)
      if (this->_lines[i][i] != 0 && i == j)
	return (i);
    return (-1);
  }

  // Calculate trace
  const TYPE	trace() const
  {
    TYPE res = 0;
    if (this->_i != this->_j)
      {
	std::cerr << "Matrix isn't nxn, cannot give trace" << std::endl;
	return (0);
      }
    for (unsigned int i = 0; i < this->_i; i++)
      res += this->_lines[i][i];
    return (res);
  }

  // Transpose
  Matrix&	transpose()
  {
    if (this->_i != this->_j)
      {
	std::cerr << "Matrix isn't nxn, cannot transpose" << std::endl;
	return (*this);
      }
    TYPE tmp = 0;
    for (unsigned int i = 0; i < this->_i; i++)
      for (unsigned int j = 0; j < i; j++)
	{
	  tmp = this->_lines[i][j];
	  this->_lines[i][j] = this->_lines[j][i];
	  this->_lines[j][i] = tmp;
	}
    return (*this);
  }

  Matrix	transpose() const
  {
    Matrix mat(this->_j, this->_i);
    for (unsigned int i = 0; i < this->_i; i++)
      for (unsigned int j = 0; j < this->_j; j++)
	mat._lines[i][j] = this->_lines[j][i];
    return (mat);
  }

  // Display matrix
  void print(std::ostream &os) const
  {
  for (unsigned int i = 0; i < this->_i; ++i)
    {
      for (unsigned int j = 0; j < this->_j; ++j)
	os << this->_lines[i][j] << " ";
      os << std::endl;
    }
  }

protected:
  typedef TYPE* _line;
  _line* _lines;
  unsigned int _i;	// Number of lines
  unsigned int _j;	// Number of column
};

template <class TYPE>
inline std::ostream& operator<<(std::ostream& ostr, const Matrix<TYPE>& stream)
{
  stream.print(ostr);
  return (ostr);
}


#endif		// MATRIX_HH