/* File: matrix.C      Copyright glh 10/21/92 */
#include  "matrix.H"

template<class T>
Matrix<T>::Matrix(int rdim, int cdim)
{
  row_dim = rdim;
  col_dim = cdim;
  element = new T* [row_dim];
  for (int i=0; i<rdim; i++)
    element[i] = new T[col_dim];
}

template<class T>
Matrix<T>::Matrix(int rdim, int cdim, const T& initval)
{
  row_dim = rdim;
  col_dim = cdim;
  element = new T* [row_dim];
  for (int i=0; i<rdim; i++)
    element[i] = new T[col_dim];
  // initialize the matrix elements
  for (i=0; i<rdim; i++)
    for (int j=0; j<cdim; j++)
      element[i][j] = initval;
}

template<class T>
Matrix<T>::Matrix(int rdim, int cdim, T* initval)
{
  row_dim = rdim;
  col_dim = cdim;
  element = new T* [row_dim];
  for (int i=0; i<rdim; i++)
    element[i] = new T[col_dim];
  // initialize the matrix elements 
  for (i=0; i<rdim; i++)
    for (int j=0; j<cdim; j++)
      element[i][j] = *(initval + (cdim * i + j));
}

template<class T>
Matrix<T>::Matrix(const Matrix<T>& m)
{
  row_dim = m.row_dim;   col_dim = m.col_dim;
  element = new T* [row_dim];
  for (int i=0; i<row_dim; i++)
    element[i] = new T[col_dim];
  for (i=0; i<row_dim; i++)
    for (int j=0; j<col_dim; j++)
      element[i][j] = m.element[i][j];
}

template<class T>
Matrix<T>::~Matrix()
{
  for (int i=0; i<row_dim; i++)
    delete [] element[i];
  delete [] element;
}

template<class T>
Matrix<T>& Matrix<T>::operator=(const Matrix<T>& m)
{
  if (element != m.element) {  // check if matrix is being assigned to itself
    for (int i=0; i<row_dim; i++)  // get rid of the old matrix
      delete [] element[i];
    delete [] element;
    row_dim = m.row_dim;   col_dim = m.col_dim;
    element = new T* [row_dim];
    for (i=0; i<row_dim; i++)
      element[i] = new T[col_dim];
    for (i=0; i<row_dim; i++)
      for (int j=0; j<col_dim; j++)
        element[i][j] = m.element[i][j];
  }
  return *this;
}

template<class T>
ostream& operator<<(ostream& os, const Matrix<T>& m) 
{
  for (int i=0; i<m.row_dim; i++)  {
    for (int j=0; j<m.col_dim; j++)  {
      os << "element(" << i << "," << j << ") = ";
      os << m.element[i][j] << endl;
    }
  }
  return os;
}
 
template<class T>
Matrix<T>& operator+(const Matrix<T>& m1, const Matrix<T>& m2)
{
  assert((m1.row_dim == m2.row_dim) && (m1.col_dim == m2.col_dim));
  Matrix<T>* result = new Matrix<T>(m1.row_dim, m1.col_dim);
  for (int i=0; i<m1.row_dim; i++)
    for (int j=0; j<m1.col_dim; j++)
      result->element[i][j] = m1.element[i][j] + m2.element[i][j];
  return *result;
}

template<class T>
Matrix<T>& operator-(const Matrix<T>& m1, const Matrix<T>& m2)
{
  assert((m1.row_dim == m2.row_dim) && (m1.col_dim == m2.col_dim));
  Matrix<T>* result = new Matrix<T>(m1.row_dim, m1.col_dim);
  for (int i=0; i<m1.row_dim; i++)
    for (int j=0; j<m1.col_dim; j++)
      result->element[i][j] = m1.element[i][j] - m2.element[i][j];
  return *result;
}

template<class T>
Matrix<T>& operator-(const Matrix<T>& m)
{
  Matrix<T>* result = new Matrix<T>(m.row_dim, m.col_dim);
  for (int i=0; i<m.row_dim; i++)
    for (int j=0; j<m.col_dim; j++)
      result->element[i][j] = -m.element[i][j]; 
  return *result;
}

template<class T>
Matrix<T>& operator*(const Matrix<T>& m1, const Matrix<T>& m2)
{
  assert(m1.col_dim == m2.row_dim);
  Matrix<T>* result = new Matrix<T>(m1.row_dim, m2.col_dim);
  for (int i=0; i<m1.row_dim; i++)
    for (int j=0; j<m2.col_dim; j++) {
      result->element[i][j] = 0;
      for (int k=0; k<m1.col_dim; k++)
        result->element[i][j] += m1.element[i][k] * m2.element[k][j];
    }
  return *result;
}

template<class T>
Matrix<T>& operator*(const T& scalar, const Matrix<T>& m)
{
  Matrix<T>* result = new Matrix<T>(m.row_dim, m.col_dim);
  for (int i=0; i<m.row_dim; i++)
    for (int j=0; j<m.col_dim; j++)
      result->element[i][j] = scalar * m.element[i][j]; 
  return *result;
}

template<class T>
Matrix<T>& Transpose(const Matrix<T>& m)
{
  Matrix<T>* result = new Matrix<T>(m.row_dim, m.col_dim);
  for (int i=0; i<m.row_dim; i++)
    for (int j=0; j<m.col_dim; j++)
      result->element[i][j] = m.element[j][i];
  return *result;
}