// matrix.h

#ifndef matrix_H
#define matrix_H

#undef  trf_info
#define trf_error
 
#include "smatrix.h"

template <class T>
class matrix : public array<T> {
private:
  int _nrow;
  int _ncol;
  int ij( int i, int j ) const { return _ncol*i+j; };
  void ostr(std::ostream& stream) const;
public:
  matrix(int M,int N) : array<T>(M*N) { _nrow = M; _ncol = N; };
  virtual ~matrix(void) { };
  matrix(const matrix<T>& mat);
  matrix(const smatrix<T>& sma);
  matrix<T> &operator=(const matrix<T>& mat);
  int nrow() const { return _nrow; };
  int ncol() const { return _ncol; };
  T& operator()( int i, int j ) const { return data[ij(i,j)]; };
  matrix<T> transpose( );
  matrix<T>& operator+=(const matrix<T>& mat);
  matrix<T>& operator-=(const matrix<T>& mat);
  matrix<T>& operator*=(const T& tval);
#if defined(DEFECT_GUIDING_DECLS)
  friend matrix<T> operator+(const matrix<T>& mat1,const matrix<T>& mat2);
  friend matrix<T> operator-(const matrix<T>& mat1,const matrix<T>& mat2);
  friend smatrix<T> operator%(const matrix<T>& mat,const smatrix<T>& sma);
#elif defined(DEFECT_NO_EXPLICIT_QUALIFICATION)
  template<class S>
  friend matrix<S> operator+(const matrix<S>& mat1,const matrix<S>& mat2);
  template<class S>
  friend matrix<S> operator-(const matrix<S>& mat1,const matrix<S>& mat2);
  template<class S>
  friend smatrix<S> operator%(const matrix<S>& mat,const smatrix<S>& sma);
#else
  friend matrix<T> operator+ <> (const matrix<T>& mat1,const matrix<T>& mat2);
  friend matrix<T> operator- <> (const matrix<T>& mat1,const matrix<T>& mat2);
  friend smatrix<T> operator%<> (const matrix<T>& mat, const smatrix<T>& sma);
#endif
};

// copier
template <class T>
inline matrix<T>::matrix(const matrix<T>& mat)
: array<T>(mat)
{
#ifdef trf_info
 std::cout << "copying a matrix\n";
#endif
 _nrow = mat._nrow;
 _ncol = mat._ncol;
}

// assignment
template <class T>
inline matrix<T> &matrix<T>::operator=(const matrix<T>& mat)
{
#ifdef trf_info
  std::cout << "assigning a matrix\n";
#endif
  assert( _ncol == mat._ncol );
  if ( _ncol == mat._ncol ) (array<T>&)*this = mat;
  else {
#ifdef trf_error
    std::cerr << "matrix::operator=: # columns does not match\n";
  }
#endif
  return *this;
}

// transpose
template <class T>
inline matrix<T> matrix<T>::transpose( ) {
 matrix<T> mat( _ncol, _nrow );
 for ( int irow=0; irow<mat.nrow(); irow++ ) {
  for ( int icol=0; icol<mat.ncol(); icol++ ) {
   mat(irow,icol) = operator()(icol,irow);
  }
 }
 return mat;
}
 

// +=
template <class T>
inline matrix<T> &matrix<T>::operator+=(const matrix<T>& mat)
{
  assert( _ncol == mat._ncol );
  if ( _ncol == mat._ncol ) (array<T>&)*this += mat;
  else {
#ifdef trf_error
    std::cerr << "matrix::operator+=: # columns does not match\n";
#endif
  }
  return *this;
}

// -=
template <class T>
inline matrix<T> &matrix<T>::operator-=(const matrix<T>& mat)
{
  assert( _ncol == mat._ncol );
  if ( _ncol == mat._ncol ) (array<T>&)*this -= mat;
  else {
#ifdef trf_error
    std::cerr << "matrix::operator-=: # columns does not match\n";
#endif
  }
  return *this;
}

// *=
template <class T>
inline matrix<T> &matrix<T>::operator*=(const T& t)
{
  (array<T>&)*this *= t;
  return *this;
}

template <class T>
inline matrix<T> operator*(const T& val,const matrix<T>& mat)
{
 matrix<T> lhs = mat;
 lhs *= val;
 return lhs;
}

template <class T>
inline matrix<T> operator*(const matrix<T>& mat,const T& val)
{
 matrix<T> lhs = mat;
 lhs *= val;
 return lhs;
}

// addition: matrix+matrix
template <class T>
inline matrix<T> operator+(const matrix<T>& mat1,const matrix<T>& mat2)
{
  matrix<T> sum = mat1;
  sum += mat2;
  return sum;
}

// subtraction: matrix+matrix
template <class T>
inline matrix<T> operator-(const matrix<T>& mat1,const matrix<T>& mat2)
{
  matrix<T> sum = mat1;
  sum -= mat2;
  return sum;
}

template <class T>
matrix<T> operator*(const matrix<T>& mat1,const matrix<T>& mat2);

template <class T>
nvector<T> operator*(const matrix<T>& mat1,const nvector<T>& vec2);

//
// mutiplication: matrix = smatrix*smatrix
//
template <class T>
matrix<T> operator*(const smatrix<T>& sma1,const smatrix<T>& sma2);

//
// mutiplication: smatrix*matrix
//
template <class T>
inline matrix<T> operator*(const smatrix<T>& sma1,const matrix<T>& mat2)
 { return matrix<T>(sma1)*mat2; }

//
// mutiplication: matrix*smatrix
//
template <class T>
inline matrix<T> operator*(const matrix<T>& mat1,const smatrix<T>& sma2)
 { return mat1*matrix<T>(sma2); }

// multiplication: M % S == M * S * MT
//
template <class T>
smatrix<T> operator%(const matrix<T>& mat, const smatrix<T>& sma);

#ifndef __xlC__
#include "matrix.c"
#endif

#endif