#ifndef MATRIX_CPP
#define MATRIX_CPP

#include <iostream>
#include "Matrix.h"

using namespace std;

template <class T>
Matrix<T>::Matrix() : m_nrow(0), m_ncol(0), m_data(NULL), m_rowPtr(NULL)
{
  cerr << "** Default constructor called." << endl;
}
  
template <class T>
Matrix<T>::Matrix(int nrow, int ncol) : m_nrow(nrow), m_ncol(ncol) 
{

  cerr << "** Non-default constructor called." << endl;

  m_data = new T[m_nrow * m_ncol]; // allocate data array
  m_rowPtr = new T*[m_nrow];       // allocate row ptr array

  for (int i = 0; i < m_nrow; i++) {
    m_rowPtr[i] = m_data + i * m_ncol;  // initialize row ptrs
    for (int j = 0; j < m_ncol; j++) {
      m_rowPtr[i][j] = T(0);            // intialize data
    }
  }
}

template <class T>
Matrix<T>::Matrix(const Matrix<T> &m) : m_nrow(m.m_nrow), m_ncol(m.m_ncol) {

  cerr << "** Copy constructor called." << endl;

  m_data = new T[m_nrow * m_ncol];   // allocate data array
  m_rowPtr = new T*[m_nrow];         // allocate row ptr array

  for (int i = 0; i < m_nrow; i++) {
    m_rowPtr[i] = m_data + i * m_ncol;  // initialize row ptrs
    for (int j = 0; j < m_ncol; j++)
      at(i, j) = m.at(i,j);             // copy m[i,j] to host[i,j]
  }
}

template <class T>
Matrix<T>::~Matrix() {

  cerr << "** Destructor called." << endl;

  if (m_data != NULL)
    delete [] m_data;      // delete data array
  if (m_rowPtr != NULL)
    delete [] m_rowPtr;    // delete row ptr array
}

template <class T>
Matrix<T>& Matrix<T>::operator=(const Matrix<T> &m) {

  cout << "** operator= called." << endl;
  
  // guard against self-assignment
  // Here m is a variable, so &m is the "address of" m

  if (this != &m) {

    // if lhs is different size from rhs, delete and 
    // re-allocate lhs arrays

    if (m_nrow != m.m_nrow || m_ncol != m.m_ncol) {

      if (m_data != NULL)
	delete [] m_data;

      if (m_rowPtr != NULL)
	delete [] m_rowPtr;

      m_ncol = m.m_ncol;
      m_nrow = m.m_nrow;
    
      m_data = new T[m_nrow * m_ncol];
      m_rowPtr = new T*[m_nrow];

      for (int i = 0; i < m_nrow; i++) {
	m_rowPtr[i] = m_data + i * m_ncol; 
      }
    }

    // copy data from rhs to lhs

    for (int i = 0; i < m_nrow; i++)
      for (int j = 0; j < m_ncol; j++)
	at(i,j) = m.at(i,j);
  }
  
  return *this;
}


template <class T>
T& Matrix<T>::at(int i, int j) const {
  return m_rowPtr[i][j];
}

template <class T>
Matrix<T> Matrix<T>::add(Matrix<T> &m) {

  if (m_nrow != m.m_nrow || m_ncol != m.m_ncol)
    throw "Matrix<T>::add() incompatible matrix sizes.";

  Matrix<T> ans(m_nrow, m_ncol);
  for (int i = 0; i < m_nrow; i++)
    for (int j = 0; j < m_ncol; j++)
      ans.at(i,j) = at(i,j) + m.at(i,j);
  return ans;
}

template <class T>
ostream& operator<<(ostream& sout, const Matrix<T> &m) {
  for (int i = 0; i < m.m_nrow; i++) {
    for (int j = 0; j < m.m_ncol; j++) {
      sout << m.at(i,j) << "\t";
    }
    sout << endl;
  }
  return sout;
}


template <class T>
T** Matrix<T>::cArray(int *nrow, int *ncol) {

  T* data = new T[m_nrow * m_ncol];  // allocate data array
  T** array = new T*[m_nrow];        // allocate row ptr array

  for (int i = 0; i < m_nrow; i++) {
    array[i] = data + i * m_ncol;
    for (int j = 0; j < m_ncol; j++)
      array[i][j] = at(i,j);
  }

  // nrow and ncol are int* variables, so *nrow and *ncol
  // dereference the variables (i.e. *nrow is the int pointed
  // to by nrow; ditto for ncol)

  *nrow = m_nrow;
  *ncol = m_ncol;
  
  return array;  // return row ptr array
}
 
#endif