#include <iostream.h>
#include <stdlib.h>
#include <math.h>
#define DEFAULT_ALLOC 2

//Defining ARRAYS_DEBUG turns on bounds checking and size checking
#define ARRAYS_DEBUG

/*
This file implements Vector and matrix for doubles and ints.
Some  arithmetic operations are defined.
An LU and LU backsubstitution is provided for
solving non-signular linear equations.
It still needs linear algebra operations defined.
*/
class matrix;
class imatrix;
class iVector;
class Vector
{
  friend class matrix;
  double *data;
  int len;

public:
  //constructors
  Vector();
  Vector(int n);
  
  //copy operator
  Vector(const Vector& original);
  
  //destructor
  ~Vector();
  
  //assignment operator
  Vector& operator =(const Vector &original);
  Vector& operator =(const Vector original);  //not really different or needed
  
  double& operator[](int i)const ;
  Vector operator +(const Vector& v);
  Vector operator -(const Vector& v);
  friend Vector operator*(double,Vector);
  friend Vector operator*(Vector,double);
  //friend Vector operator*(Vector&);
  
  //printing for debugging purposes
  void rprint() const;
  void print()const ;

  void resize(int n);

  friend double dot(const Vector& v1, const Vector& v2);
  friend double norm(const Vector& v);
  double norm(void)const;
  void zero(void);
  int length() const{ return len;}
  int operator==(const Vector& v)const;
  double max(void)const;
  double min(void)const;
  friend ostream& operator<<(ostream& s,Vector v);
  double mean()const;
};

class iVector
{
  friend class imatrix;
  int *data;
  int len;

public:
	//constructors
  iVector();
  iVector(int n);
  //copy operator
  iVector(const iVector& original);
  //destructor
  ~iVector();
  //assignment operator
  iVector& operator =(const iVector &original);
  iVector& operator =(const iVector original);
  int& operator[](int i)const;
  iVector operator +(const iVector& v);
  iVector operator -(const iVector& v);
  friend iVector operator*(int,iVector);
  friend iVector operator*(iVector,int);
  //friend iVector operator*(iVector&);
  void rprint() const;
  void print()const ;
  friend int dot(const iVector& v1, const iVector& v2);
  friend double norm(const iVector& v);
  int max(void)const;
  int min(void)const;
  int length() const{ return len;}
  void zero(void);
  int operator==(const iVector& v)const;
  friend ostream& operator<<(ostream& s,iVector v);
};
/*
Vector class definitions
*/
double Vector::max(void)const
{
	double max=data[0];
	for(int i=1;i<len;i++) if(data[i] > max) max=data[i];
	return max;
}
double Vector::min(void)const
{
	double min=data[0];
	for(int i=1;i<len;i++) if(data[i] < min) min=data[i];
	return min;
}
ostream& operator<<(ostream& s,Vector v)
{
	s << "(";
	for(int i=0;i<v.len-1;i++) s << v.data[i] << ", ";
	s << v.data[v.len-1]<<")"<<endl;
	return s;
}
int Vector::operator==(const Vector& v)const
{
	if(len != v.len) return 0;
	for(int i=0;i<len;i++) if(data[i]!=v.data[i]) return 0;
	return 1;
}
double& Vector::operator[](int i) const
{
#ifdef ARRAYS_DEBUG
	if( i<0 || i >= len)
	{
		cerr << "ERROR Vector::operator[](int i) const with i = " << i << " out of bounds" << endl;
		cerr << "The length of the vector is " << len << endl;
		print();
		exit(2);
	}
#endif
	return data[i];
}
Vector::~Vector(){delete [] data;}
Vector::Vector(int n)
{
	data = new double[len=n];
	if(!data)
	{
		cerr << "VECTOR Error: Vector::Vector(int n) can't allocate memory with n = " << n << endl;
		exit(2);
	}
}
Vector::Vector()
{
	data = new double[DEFAULT_ALLOC];
	len=DEFAULT_ALLOC;
	if(!data)
	{
		cerr << "VECTOR Error: Vector::Vector() can't allocate memory" << endl;
		exit(2);
	}
}
Vector operator*(double x,Vector v)
{
	int i;
	Vector ans(v.len);
	for(i=0;i<v.len;i++) ans.data[i]= x*v[i];
	return ans;
}
Vector operator*(Vector v, double x)
{
	int i;
	Vector ans(v.len);
	for(i=0;i<v.len;i++) ans.data[i]= x*v[i];
	return ans;
}
Vector Vector::operator +(const Vector& v)
{
	int i;
#ifdef ARRAYS_DEBUG
	if( len != v.len)
	{
		cerr << "VECTOR Error: You're trying to add vectors of different sizes\n";
		rprint();
		v.rprint();
		exit(2);
	}
#endif
	Vector sum(len);
	for(i=0;i<len;i++) sum.data[i] = data[i]+v.data[i];
	return sum;
}
Vector Vector::operator -(const Vector& v)
{
	int i;
#ifdef ARRAYS_DEBUG
	if( len != v.len)
	{
		cerr << "VECTOR Error: You're trying to subtract vectors of different sizes\n";
		rprint();
		v.rprint();
		exit(2);
	}
#endif
	Vector sum(len);
	for(i=0;i<len;i++) sum.data[i] = data[i]-v.data[i];
	return sum;
}
double Vector::mean()const
{
	double sum=0;
	for(int i=0;i<len;i++) sum += data[i];
	return sum/len;
}
double dot(const Vector& v1,const Vector& v2)
{
  int i;
  double s=0.0;
#ifdef ARRAYS_DEBUG
  if( v1.len != v2.len)
	 {
		cerr << "Error: Different sized vectors in dot()--- returning 0.0\n";
		v1.rprint();
		v2.rprint();
		return s;
	 }
#endif
  for(i=0;i<v1.len;i++) s += v1[i]*v2[i];
  return s;
}


void  Vector::rprint()const  //print entries on a single line
{
  int i;
  cout << "VECTOR: ";
  cout << "(";
  for(i=0;i<len-1;i++) cout << data[i] << ",";
  cout << data[len-1] << ")" << endl;
  return;
}

void Vector::print() const
{
  int i;
  cout << "VECTOR: Length = " << len << "\n";
  for(i=0;i<len;i++) cout << data[i] <<"\n";
  return;
}


Vector::Vector(const Vector& original)
{
  data=new double[len=original.len];
  int i;
  for(i=0;i<len;i++) data[i]=original.data[i];
}

void Vector::resize(int n)
{
	delete[]data;
	data = new double[len=n];
}
Vector& Vector::operator=(const Vector &original)
{
  if(this != &original)
  {
	delete [] data;
	data= new double[len=original.len];
	int i;
	for(i=0;i<len;i++) data[i]=original.data[i];
	}
  return *this;
}

Vector& Vector::operator=(const Vector original)
{
  if(this != &original)
  {
	delete [] data;
	data= new double[len=original.len];
	int i;
	for(i=0;i<len;i++) data[i]=original.data[i];
	}
  return *this;
}

double Vector::norm(void) const
{
	double s=0.0;
	for(int i=0;i<len;i++) s+= data[i]*data[i];
	return sqrt(s);
}
//I don't imagine one ever needs this form of norm()
double norm(const Vector& v)
{
	double s=0.0;
	for(int i=0; i< v.len; i++) s+= v.data[i]*v.data[i];
	return sqrt(s);
}
void Vector::zero(void)
{
	for(int i=0;i<len;i++) data[i]=0.0;
	return;
}

/**************************
iVector   class definitions
***************************/
void iVector::zero()
{
	for(int i=0;i<len;i++) data[i]=0;
	return;
}
ostream& operator<<(ostream& s,iVector v)
{
	s << "(";
	for(int i=0;i<v.len-1;i++) s << v.data[i] << ", ";
	s << v.data[v.len-1]<<")"<<endl;
	return s;
}
int iVector::max()const
{
	int max=data[0];
	for(int i=1;i<len;i++) if(data[i] > max) max=data[i];
	return max;
}
int iVector::min()const
{
	int min=data[0];
	for(int i=1;i<len;i++) if(data[i] < min) min=data[i];
	return min;
}
int iVector::operator==(const iVector& v)const
{
	if(len != v.len) return 0;
	for(int i=0;i<len;i++) if(data[i]!=v.data[i]) return 0;
	return 1;
}
int& iVector::operator[](int i) const
{
#ifdef ARRAYS_DEBUG
	if( i<0 || i >= len)
	{
		cerr << "ERROR iVector::operator[](int i) const with i = " << i << " out of bounds" << endl;
		print();
		exit(2);
	}
#endif
	return data[i];
}
iVector::~iVector(){delete [] data;}
iVector::iVector(int n)
{
	data = new int[len=n];
	if(!data)
	{
		cerr << "VECTOR Error: iVector::iVector(int n) can't allocate memory with n = " << n << endl;
		exit(2);
	}
}
iVector::iVector()
{
	data = new int[DEFAULT_ALLOC];
	len=DEFAULT_ALLOC;
	if(!data)
	{
		cerr << "VECTOR Error: iVector::iVector() can't allocate memory" << endl;
		exit(2);
	}
}
iVector operator*(int x,iVector v)
{
	int i;
	iVector ans(v.len);
	for(i=0;i<v.len;i++) ans.data[i]= x*v[i];
	return ans;
}
iVector operator*(iVector v, int x)
{
	int i;
	iVector ans(v.len);
	for(i=0;i<v.len;i++) ans.data[i]= x*v[i];
	return ans;
}
iVector iVector::operator +(const iVector& v)
{
	int i;
#ifdef ARRAYS_DEBUG
	if( len != v.len)
	{
		cerr << "VECTOR Error: You're trying to add vectors of different sizes\n";
		rprint();
		v.rprint();
		exit(2);
	}
#endif
	iVector sum(len);
	for(i=0;i<len;i++) sum.data[i] = data[i]+v.data[i];
	return sum;
}
iVector iVector::operator -(const iVector& v)
{
	int i;
#ifdef ARRAYS_DEBUG
	if( len != v.len)
	{
		cerr << "VECTOR Error: You're trying to subtract vectors of different sizes\n";
		rprint();
		v.rprint();
		exit(2);
	}
#endif
	iVector sum(len);
	for(i=0;i<len;i++) sum.data[i] = data[i]-v.data[i];
	return sum;
}

int dot(const iVector& v1,const iVector& v2)
{
  int i;
  int s=0;
#ifdef ARRAYS_DEBUG
  if( v1.len != v2.len)
	 {
		cerr << "Error: Different sized vectors in dot()--- returning 0\n";
		v1.rprint();
		v2.rprint();
		return s;
	 }
#endif
  for(i=0;i<v1.len;i++) s += v1[i]*v2[i];
  return s;
}


void  iVector::rprint()const  //print entries on a single line
{
  int i;
  cout << "VECTOR: ";
  cout << "(";
  for(i=0;i<len-1;i++) cout << data[i] << ",";
  cout << data[len-1] << ")" << endl;
  return;
}

void iVector::print() const
{
  int i;
  cout << "VECTOR: Length = " << len << "\n";
  for(i=0;i<len;i++) cout << data[i] <<"\n";
  return;
}


iVector::iVector(const iVector& original)
{
  data=new int[len=original.len];
  int i;
  for(i=0;i<len;i++) data[i]=original.data[i];
}

iVector& iVector::operator=(const iVector &original)
{
  if(this != &original)
  {
	delete [] data;
	data= new int[len=original.len];
	int i;
	for(i=0;i<len;i++) data[i]=original.data[i];
	}
  return *this;
}

iVector& iVector::operator=(const iVector original)
{
  if(this != &original)
  {
	delete [] data;
	data= new int[len=original.len];
	int i;
	for(i=0;i<len;i++) data[i]=original.data[i];
	}
  return *this;
}
double norm(const iVector& v)
{
	int s=0;
	for(int i=0; i< v.len; i++) s+= v.data[i]*v.data[i];
	return sqrt((double)s);
}


/**********************
MATRIX
**********************/

class matrix
{
	Vector* m;
	int rows,cols;
	public:
	//constructor
	matrix();
	matrix(const int r,const int c);
	//copy
	matrix(const matrix &s);
	//destructor
	~matrix();
	//assignment operator
	matrix& operator =(const matrix& s);
	void print() const;
	void pprint() const;
	//select a row
	Vector& operator[](const int i);
	//arithmetic operations
	Vector operator*(const Vector&);
	friend matrix  operator*(const double, const matrix);
	friend matrix operator*(const matrix, const double);
	matrix operator*(const matrix a);
	matrix operator+(const matrix& a);
	matrix operator-(const matrix& a);
	friend void ludcmp(matrix& a,iVector& indx,double &d);
	friend void lubksb(matrix&a,iVector& indx,Vector&b);
	double rowsum(int i) const;
	void zero(void);
	void identity(void);
	friend ostream& operator<<(ostream& s,matrix m);
};

class imatrix
{
	iVector* m;
	int rows,cols;
	public:
	//constructor
	imatrix();
	imatrix(const int r,const int c);
	//copy
	imatrix(const imatrix &s);
	//destructor
	~imatrix();
	//assignment operator
	imatrix& operator =(const imatrix& s);
	void print() const;
	//select a row
	iVector& operator[](const int i);
	//arithmetic operations
	iVector operator*(const iVector&);
	friend imatrix  operator*(const int, const imatrix);
	friend imatrix operator*(const imatrix, const int);
	imatrix operator+(const imatrix& a);
	imatrix operator-(const imatrix& a);
	int rowsum(int i)const ;
	void zero(void);
	friend ostream& operator<<(ostream& s,imatrix m);

};

/*
matrix class definitions
*/
matrix::matrix()
{
	m = new Vector[DEFAULT_ALLOC];
	if(!m)
	{
		cerr << "MATRIX Error: matrix::matrix() can't allocate memory" << endl;
		exit(2);
	}
	rows=cols=DEFAULT_ALLOC;
	for(int i=0;i<rows;i++)
	{
		Vector v;
		m[i]= v;
	}
}
matrix::matrix(const int r,const int c)
{
	m= new Vector[r];
	if(!m)
	{
		cerr << "MATRIX Error: matrix::matrix(int r, int c) can't allocate memory" << endl;
		cerr << "r = " << r << "c = " << c << endl;
		exit(2);
	}
	rows=r;
	cols=c;
	for(int i=0;i<r;i++)
	{
		Vector v(cols);
		m[i]=v;
	}
}

matrix::matrix(const matrix &s)
{
	int i;
	rows=s.rows;
	m = new Vector[rows];
	if(!m)
	{
		cerr << "MATRIX Error: matrix::matrix(const matrix &s) can't allocate memory" << endl;
		s.print();
		exit(2);
	}
	cols =s.cols;
	for(i=0;i<rows;i++)
	{
	  m[i]=s.m[i];
	}
}

matrix::~matrix()
{
	delete [] m;
}
matrix& matrix::operator =(const matrix &s)
{
	if(this != &s)
	{
		delete []m;
		rows= s.rows;
		cols=s.cols;
		m = new Vector[rows];
		for(int i=0;i<rows;i++) m[i]=s.m[i];
	}
	return *this;
}
void matrix::print()const
{
	cout << "MATRIX: Rows = " << rows << " Cols = " << cols << endl;
	for(int i=0;i<rows;i++)
	{
		 m[i].rprint();
	}
}
void matrix::pprint()const
{
	for (int i=0;i<rows;i++)
	{
		for(int j=0;j<cols;j++) cout << m[i][j] << " ";
		cout << endl;
	}
	return;
}
ostream& operator<<(ostream& s,matrix m)
{
	for(int i=0; i<m.rows;i++) s << m[i];
	return s;
}
Vector& matrix::operator[](const int i)
{
#ifdef ARRAYS_DEBUG
	if(i<0 || i>=rows)
	{
		cerr << "MATRIX Error: select row operator index " << i << " out of range\n";
		print();
		exit(2);
	}
#endif
	return m[i];
}
Vector matrix::operator*(const Vector& v)
{
	int i,j;
#ifdef ARRAYS_DEBUG
	if(cols != v.len)
	{
		cerr << "MATRIX Error: multiplying matrix times Vector with wrong sizes\n";
		print();
		v.rprint();
		exit(2);
	}
#endif
	Vector ans(rows);
	for(i=0;i<rows;i++)
	{
		ans.data[i]=0.0;
		for(j=0;j<cols;j++) ans.data[i] += m[i][j]*v.data[j];
	}
	return ans;
}
matrix operator*(const double x,const matrix s)
{
	matrix ans(s.rows,s.cols);
	for(int i=0;i<ans.rows;i++)
	  {
		ans.m[i]= x*s.m[i];
	  }
	return ans;
}
matrix operator*(const matrix s,const double x)
{
	int i;
	matrix ans(s.rows,s.cols);
	for(i=0;i<ans.rows;i++)
	  {
		 ans.m[i]=x*s.m[i];
	  }
	return ans;
}
matrix matrix::operator*(const matrix a)
{
#ifdef ARRAYS_DEBUG
	if(cols != a.rows)
	{
		cerr << "MATRIX Error: matrix matrix::operator*(const matrix a):" << endl;
		cerr << "matricies are the wrong size " << cols << " " << a.rows << endl;
		print();
		a.print();
		exit(2);
	}
#endif
	matrix ans(rows,a.cols);
	for(int i=0;i<rows;i++)
	{
		for(int j=0;j<a.cols;j++)
		{
			ans.m[i][j]=0.0;
			for(int k=0;k<cols;k++)
			{
				ans.m[i][j] += m[i][k]*a.m[k][j];
			}
		}
	}
	return ans;
}
matrix matrix::operator+(const matrix& a)
{
	int i,j;
#ifdef ARRAYS_DEBUG
	if(rows!= a.rows || cols != a.cols)
	{
		cerr << "MATRIX Error: you're trying to add matricies of different sizes\n";
		print();
		a.print();
		exit(2);
	}
#endif
	matrix ans(a.rows,a.cols);
	for(i=0;i<a.rows;i++)
	{
		for(j=0;j<a.cols;j++)
		  {
			ans.m[i][j] = m[i][j] + a.m[i][j];  //faster than assigning vectors?
		}
	}
	return ans;
}
matrix matrix::operator-(const matrix& a)
{
	int i,j;
#ifdef ARRAYS_DEBUG
	if(rows != a.rows || cols != a.cols)
	{
		cerr << "ERROR: you're trying to subtract matricies of different sizes\n";
		print();
		a.print();
		exit(2);
	}
#endif
	matrix ans(rows,cols);
	for(i=0;i<rows;i++)
	{
		for(j=0;j<cols;j++)
		ans.m[i][j] = m[i][j] - a.m[i][j];
	}
	return ans;
}

double matrix::rowsum(int i) const
{
#ifdef ARRAYS_DEBUG
	if(i<0 || i >= rows)
	{
		cerr << "ERROR: matrix::rowsum(int i) row number " << i << "out of bounds" << endl;
		print();
		exit(2);
	}
#endif
	double sum=0;
	for(int k=0;k<cols;k++) sum += m[i][k];
	return sum;
}
void matrix::zero(void)
{
	for(int i=0;i<rows;i++)
	{
		for(int j=0;j<cols;j++)
		{
			m[i][j]=0.0;
		}
	}
	return;

}

void matrix::identity(void)
{
#ifdef ARRAYS_DEBUG
	if(rows!=cols)
	{
		cerr << "MATRIX Error: matrix::identity(): Matrix not square\n";
		print();
		exit(2);
	}
#endif
	zero();
	for(int i=0;i<rows;i++) m[i][i]=1.0;
	return;
}

#define TINY 1.0e-20;
void ludcmp(matrix& a, iVector& indx,double& d)
{
	int i,imax,j,k;
	double  big,dum,sum,temp;
	int n=a.rows;
	Vector vv(n);
	if(a.rows != a.cols)
	{
		cerr << "MATRIX ERROR: ludcmp(). Matrix isn't square." << endl;
		a.print();
		exit(2);
	}
	d=1.0;
	for (i=0;i<n;i++)
	{
		big=0.0;
		for (j=0;j<n;j++) if ((temp=fabs(a[i][j])) > big) big=temp;
		if (big == 0.0) cerr << "Singular matrix in routine LUDCMP" << endl;
		vv[i]=1.0/big;
	}
	for (j=0;j<n;j++)
	{
		for (i=0;i<j;i++)
		{
			sum=a[i][j];
			for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
		}
		big=0.0;
		for (i=j;i<n;i++)
		{
			sum=a[i][j];
			for (k=0;k<j;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
			if ( (dum=vv[i]*fabs(sum)) >= big)
			{
				big=dum;
				imax=i;
			}
		}
		if (j != imax)
		{
			for (k=0;k<n;k++)
			{
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			d = -(d);
			vv[imax]=vv[j];
		}
		indx[j]=imax;
		if (a[j][j] == 0.0) a[j][j]=TINY;
		if (j != n-1) {
			dum=1.0/(a[j][j]);
			for (i=j+1;i<n;i++) a[i][j] *= dum;
		}
	}
}
#undef TINY
void lubksb(matrix& a,iVector& indx,Vector& b)
{
	int i,ip,j;
	double sum;
	int n=a.rows;
	for (i=0;i<n;i++)
	{
		ip=indx[i];
		sum=b[ip];
		b[ip]=b[i];
		for (j=0;j<=i-1;j++) sum -= a[i][j]*b[j];
		b[i]=sum;
	}
	for (i=n-1;i>=0;i--)
	{
		sum=b[i];
		for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
		b[i]=sum/a[i][i];
	}
}

/***********************
imatrix   class definitions
************************/
imatrix::imatrix()
{
	m = new iVector[DEFAULT_ALLOC];
	if(!m)
	{
		cerr << "MATRIX Error: imatrix::imatrix() can't allocate memory" << endl;
		exit(2);
	}
	rows=cols=DEFAULT_ALLOC;
	for(int i=0;i<rows;i++)
	{
		iVector v;
		m[i]= v;
	}
}
imatrix::imatrix(const int r,const int c)
{
	m= new iVector[r];
	if(!m)
	{
		cerr << "MATRIX Error: imatrix::imatrix(int r, int c) can't allocate memory" << endl;
		cerr << "r = " << r << "c = " << c << endl;
		exit(2);
	}
	rows=r;
	cols=c;
	for(int i=0;i<r;i++)
	{
		iVector v(cols);
		m[i]=v;
	}
}

imatrix::imatrix(const imatrix &s)
{
	int i;
	rows=s.rows;
	m = new iVector[rows];
	if(!m)
	{
		cerr << "MATRIX Error: imatrix::imatrix(const imatrix &s) can't allocate memory" << endl;
		s.print();
		exit(2);
	}
	cols =s.cols;
	for(i=0;i<rows;i++)
	{
	  m[i]=s.m[i];
	}
}

imatrix::~imatrix()
{
	delete [] m;
}
imatrix& imatrix::operator =(const imatrix &s)
{
	if(this != &s)
	{
		delete []m;
		rows= s.rows;
		cols=s.cols;
		m = new iVector[rows];
		for(int i=0;i<rows;i++) m[i]=s.m[i];
	}
	return *this;
}
void imatrix::print()const
{
	cout << "MATRIX: Rows = " << rows << "Cols = " << cols << endl;
	for(int i=0;i<rows;i++)
	{
		 m[i].rprint();
	}
}
iVector& imatrix::operator[](const int i)
{
#ifdef ARRAYS_DEBUG
	if(i<0 || i>=rows)
	{
		cerr << "MATRIX Error: select row operator index " << i << " out of range\n";
		print();
		exit(2);
	}
#endif
	return m[i];
}
iVector imatrix::operator*(const iVector& v)
{
	int i,j;
#ifdef ARRAYS_DEBUG
	if(cols != v.len)
	{
		cerr << "MATRIX Error: multiplying imatrix times iVector with wrong sizes\n";
		print();
		v.rprint();
		exit(2);
	}
#endif
	iVector ans(rows);
	for(i=0;i<rows;i++)
	{
		ans.data[i]=0;
		for(j=0;j<cols;j++) ans.data[i] += m[i][j]*v.data[j];
	}
	return ans;
}
imatrix operator*(const int x,const imatrix s)
{
	imatrix ans(s.rows,s.cols);
	for(int i=0;i<ans.rows;i++)
	  {
		ans.m[i]= x*s.m[i];
	  }
	return ans;
}
imatrix operator*(const imatrix s,const int x)
{
	int i;
	imatrix ans(s.rows,s.cols);
	for(i=0;i<ans.rows;i++)
	  {
		 ans.m[i]=x*s.m[i];
	  }
	return ans;
}

imatrix imatrix::operator+(const imatrix& a)
{
	int i,j;
#ifdef ARRAYS_DEBUG
	if(rows!= a.rows || cols != a.cols)
	{
		cerr << "MATRIX Error: you're trying to add matricies of different sizes\n";
		print();
		a.print();
		exit(2);
	}
#endif
	imatrix ans(a.rows,a.cols);
	for(i=0;i<a.rows;i++)
	{
		for(j=0;j<a.cols;j++)
		  {
			ans.m[i][j] = m[i][j] + a.m[i][j];  //faster than assigning vectors?
		}
	}
	return ans;
}

imatrix imatrix::operator-(const imatrix& a)
{
	int i,j;
#ifdef ARRAYS_DEBUG
	if(rows != a.rows || cols != a.cols)
	{
		cerr << "ERROR: you're trying to subtract matricies of different sizes\n";
		print();
		a.print();
		exit(2);
	}
#endif
	imatrix ans(rows,cols);
	for(i=0;i<rows;i++)
	{
		for(j=0;j<cols;j++)
		ans.m[i][j] = m[i][j] - a.m[i][j];
	}
	return ans;
}

int imatrix::rowsum(int i) const
{
#ifdef ARRAYS_DEBUG
	if(i<0 || i >= rows)
	{
		cerr << "iMatrix ERROR: row number " << i << "out of bounds" << endl;
		print();
		exit(2);
	}
#endif
	int sum=0;
	for(int k=0;k<cols;k++) sum += m[i][k];
	return sum;
}

void imatrix::zero(void)
{
	for(int i=0;i<rows;i++)
	{
		for(int j=0;j<cols;j++)
		{
			m[i][j]=0;
		}
	}
	return;

}
ostream& operator<<(ostream& s,imatrix m)
{
	for(int i=0; i<m.rows;i++) s << m[i];
	return s;
}