#ifndef __P4840_MATRIX_H__
#define __P4840_MATRIX_H__

#include <iostream>

namespace p4840 
{

template<class T>
class Matrix
{
  private:
    T* _data;
    size_t num_rows, num_cols;

    void init()
      {
        _data = new T[num_rows*num_cols];
      }

  public:

//  Constructors
    Matrix(): num_rows(1), num_cols(1) 
      { init(); } 

#ifdef USE_FORTRAN_STYLE
    Matrix(size_t ncols, size_t nrows): num_rows(nrows),num_cols(ncols) 
      { init(); }
#else
    Matrix(size_t nrows, size_t ncols): num_rows(nrows),num_cols(ncols) 
      { init(); }
#endif

    Matrix(const Matrix<T>& m): num_rows(m.rows()),num_cols(m.cols())
      {
        init();
        for(int i = 0; i < num_rows*num_cols; i++) _data[i] = m[i];
      }

//  Operators

#ifdef USE_FORTRAN_STYLE
    T operator()(size_t col, size_t row) const; // Fortran style
      { return _data[row*num_cols+col]; }

    T& operator()(size_t col, size_t row);
      { return _data[row*num_cols+col]; }
#else
    T operator()(size_t row, size_t col) const // C style
      { return _data[row*num_cols+col]; }

    T& operator()(size_t row, size_t col)
      { return _data[row*num_cols+col]; }
#endif
    T operator[](size_t index) const
      { return _data[index]; }
    T& operator[](size_t index)
      { return _data[index]; }

//  Asignment
//  Missing substraction, multiplication and divition by constant...

    Matrix& operator=(const Matrix<T> &m)
      {
        delete[](_data); 
        num_cols = m.num_rows; 
        num_rows = m.num_cols; 
        init(); 
        for(int i = 0; i < num_rows*num_cols; i++) _data[i] = m[i];
        return *this;
      }
        
    Matrix& operator=(const T& v)
      {
        for(int i = 0; i < num_rows*num_cols; i++) _data[i] = v;
        return *this;
      }

    Matrix& operator+=(const Matrix<T> &m)
      {
        if(m.rows() != rows() || m.cols() != cols()) 
        { 
           std::cerr << "Matrix dimensions are different\n"; 
           return *this;
        }
        for(int i = 0; i < num_rows*num_cols; i++) _data[i] += m[i];
        return *this;
      }
          
    Matrix& operator*=(T v)
      {
        for(int i = 0; i < num_rows*num_cols; i++) _data[i] *= v;
        return *this;
      }
          
    Matrix operator+(const Matrix<T> &m)
      {
        Matrix<T> ret_val;

        if(m.rows() != rows() || m.cols() != cols()) 
        { 
           std::cerr << "Matrix dimensions are different\n"; 
           return *this;
        }
        ret_val = *this;
        ret_val += m;
        return ret_val;
      }
          
    Matrix operator*(T v)
      {
        Matrix<T> ret_val;
        for(int i = 0; i < num_rows*num_cols; i++) ret_val[i] = _data[i]*v;
        return ret_val;
      }

//  Methods
//  Missing: invert, rotate, randomize, trace, determinant,...

    int rows() const { return num_rows; }
    int cols() const { return num_cols; }

//  Memory management
#ifdef USE_FORTRAN_STYLE
    Matrix<T>& resize(size_t ncols, size_t nrows)
      { 
        delete[](_data);  
        num_cols = ncols; 
        num_rows = nrows; 
        init(); 
        return *this;  
      }
#else
    Matrix<T>& resize(size_t nrows, size_t ncols)
      { 
        delete[](_data);  
        num_cols = ncols; 
        num_rows = nrows; 
        init(); 
        return *this;  
      }
#endif

//  Destructor
    ~Matrix<T>() { delete[](_data); }
};


} // namespace p4840

#endif // __P4840_MATRIX_H