#ifndef __MATRIX_H__
#define __MATRIX_H__

/* $Id: Matrix.H,v 1.4 91/04/16 18:03:17 bjaspan Exp $ */

#include <assert.h>  /* for abort() */
#include <stream.h>

class Matrix {
public:
     /* Constructors */
     Matrix();
     Matrix(int, int);
     Matrix(int, int, int *);
     Matrix(int, int, double *);
     Matrix(int, int, int ...);
     Matrix(int, int, double ...);
     Matrix(const Matrix&);
     Matrix(const Matrix&,int,int,int,int);
     ~Matrix();

     /* Member functions */
     Matrix T() const;
     inline int Rows() const;
     inline int Cols() const;
     inline int SameSize(const Matrix&) const;
     inline double& operator()(int, int) const;
     inline double *Row(int) const;
     void Print() const;

     /* Overloaded operators */
     friend Matrix operator+(const Matrix&, const Matrix&);
     friend Matrix operator-(const Matrix&, const Matrix&);
     friend Matrix operator*(const Matrix&, const Matrix&);
     friend Matrix operator*(double, const Matrix&);
     Matrix& operator+=(const Matrix&);
     Matrix& operator-=(const Matrix&);
     inline friend Matrix& operator*=(Matrix&, const Matrix&);
     inline friend Matrix operator/(const Matrix&, double);
     friend int operator==(const Matrix&, const Matrix&);
     friend Matrix operator|(const Matrix&, const Matrix&);
     friend Matrix operator&(const Matrix&, const Matrix&);
     Matrix& operator=(const Matrix&);
     friend ostream& operator<<(ostream&, const Matrix&);

     /* Special matrix creators */
     static inline const Matrix& NAM();

     /* Matrix mutators */
     void Identify();

private:
     void Init(const int, const int);

     static Matrix *M_NAM;
     static Matrix *MakeNAM();
     
     int M_m, M_n, nam;
     double *d;
};

inline int Matrix::Rows() const
{
     return M_m;
}

inline int Matrix::Cols() const
{
     return M_n;
}

inline Matrix& operator*=(Matrix& m1, const Matrix& m)
{
   return(m1 = m1*m);
}

inline Matrix operator/(const Matrix& m, double k)
{
   return((1.0/k)*m);
}

inline double& Matrix::operator()(int m, int n) const
{
     /* My kingdom for exception handling ... */
     if (m > M_m || n > M_n || m < 1 || n < 1)
	  abort();
     
     return d[M_n*(m-1) + n - 1];
}

inline double *Matrix::Row(int r) const
{
   return(d+(r-1)*M_n);
}

inline const Matrix& Matrix::NAM()
{
     return *Matrix::M_NAM;
}

inline int Matrix::SameSize(const Matrix& m) const
{
     return (M_m == m.M_m && M_n == m.M_n);
}

/* DO NOT ADD ANYTHING AFTER THIS #endif */
#endif /* __MATRIX_H__ */