// MatrixD.H
// Given up using a Matrix<T> approach for now: too many 
// aggravating difficulties.
// News group response suggests looking at the "STL"
// I ftp'd the STL from ftp.cs.rpi.edu today (to lake-o).
// Will look at it later.

// See also the analogous int class in MatrixI.H

//	JR
//	11/3/95  Tested constructors.
//		 Added member "rows()".
//		       operators +=,+,-=,-.
//		       const operator() (int,int) const
//	11/4/95  Qi tested. Found error in copy constructors.
//		 Reversed order of deletion in destructor.

#ifndef MATRIXD_H
#define MATRIXD_H
#include <assert.h>
#include <math.h>

#include "MatrixI.H"


class MatrixD
{
   private:
	int nrows;
	int ncols;
	double **mat;
   public:
	MatrixD()
	: nrows(0), ncols(0), mat(0)
	{}
	MatrixD(int nr,int nc);
	MatrixD(int nr,int nc,double val);
	MatrixD(int n);
	~MatrixD();
	MatrixD(const MatrixD& m);
	MatrixD(const MatrixI& m);
	MatrixD& operator=(const MatrixD& m)
	{	if(&m==this) return *this;
		if(nrows!=m.nrows || ncols!=m.ncols)
		{	for(int i=0;i<nrows;i++) delete[] mat[i];
			delete[] mat;
			nrows=m.nrows;
			ncols=m.ncols;
			mat = new double*[nrows];
			for(i=0;i<nrows;i++) 
			{	mat[i] = new double[nrows];
				for(int j=0;j<ncols;j++) mat[i][j]=m.mat[i][j];
			}
		}
		else
		{	for(int i=0;i<nrows;i++) 
				for(int j=0;j<ncols;j++) 
					mat[i][j]=m.mat[i][j];
		}
		return *this;
	}
	MatrixD& operator*=(const double& d)
	{	for(int i=0;i<nrows;i++) 
			for(int j=0;j<ncols;j++) 
				mat[i][j]*=d;
		return *this;
	}
	MatrixD& operator/=(const double& d)
	{	for(int i=0;i<nrows;i++) 
			for(int j=0;j<ncols;j++) 
				mat[i][j]/=d;
		return *this;
	}

	double& operator () (int i,int j)
	{	assert( 0<=i && i<nrows && 0<=j && j<ncols );
		return mat[i][j];
	}
	const double& operator () (int i,int j) const
	{	assert( 0<=i && i<nrows && 0<=j && j<ncols );
		return mat[i][j];
	}
	double largestMag() const
	{	double largest=0.0;
		for(int i=0;i<nrows;i++) 
			for(int j=0;j<ncols;j++) 
				if(fabs(mat[i][j])>largest)
					largest=fabs(mat[i][j]);
		return largest;
	}
	int rows() const
	{ 	return nrows;
	}			

	MatrixD& operator+=(const MatrixD& m);
	MatrixD& operator-=(const MatrixD& m);

	friend MatrixD operator+(const MatrixD& a,
				 const MatrixD& b);
	friend MatrixD operator-(const MatrixD& a,
				 const MatrixD& b);
			
	friend MatrixD operator*(const MatrixD& a,
				 const MatrixD& b);
	friend ostream& operator << (ostream& out, const MatrixD& m);

};

MatrixD::MatrixD(int nr,int nc)
: nrows(nr), ncols(nc)
{	mat = new double*[nrows];
	for(int i=0;i<nrows;i++) 
	{	mat[i] = new double[ncols];
		for(int j=0;j<ncols;j++) mat[i][j]=0.0;
	}

}

MatrixD::MatrixD(int nr,int nc,double val)
: nrows(nr), ncols(nc)
{	mat = new double*[nrows];
	for(int i=0;i<nrows;i++) 
	{	mat[i] = new double[ncols];
		for(int j=0;j<ncols;j++) mat[i][j]=val;
	}
}

MatrixD::MatrixD(int n)
: nrows(n), ncols(n)
{	mat = new double*[n];
	for(int i=0;i<n;i++) 
	{	mat[i] = new double[n];
		for(int j=0;j<n;j++) mat[i][j]=0.0;
	}
}

MatrixD::MatrixD(const MatrixD& m)
: nrows(m.nrows), ncols(m.ncols)
{	mat = new double*[nrows];
	for(int i=0;i<nrows;i++) 
	{	mat[i] = new double[ncols];
		for(int j=0;j<ncols;j++) mat[i][j]=m.mat[i][j];
	}
}

MatrixD::MatrixD(const MatrixI& m)
: nrows(m.nrows), ncols(m.ncols)
{	mat = new double*[nrows];
	for(int i=0;i<nrows;i++) 
	{	mat[i] = new double[ncols];
		for(int j=0;j<ncols;j++) mat[i][j]=m.mat[i][j];
	}
}

MatrixD::~MatrixD()
{	for(int i=nrows-1;i>=0;i--) delete[] mat[i];
	delete[] mat;
}

MatrixD& MatrixD::operator+=(const MatrixD& m)
{
	assert( (nrows==m.nrows) && (ncols==m.ncols) );
	for(int i=0;i<nrows;i++)
		for(int j=0;j<ncols;j++) (*this)(i,j)+=m(i,j);
	return *this;
}

MatrixD operator+(const MatrixD& m1,const MatrixD& m2)
{
	MatrixD tmp(m1);
	return (tmp+=m2);
}

MatrixD& MatrixD::operator-=(const MatrixD& m)
{
	assert( (nrows==m.nrows) && (ncols==m.ncols) );
	for(int i=0;i<nrows;i++)
		for(int j=0;j<ncols;j++) (*this)(i,j)-=m(i,j);
	return *this;
}

MatrixD operator-(const MatrixD& m1,const MatrixD& m2)
{
	MatrixD tmp(m1);
	return (tmp-=m2);
}

MatrixD operator*(const MatrixD& a,
	 	  const MatrixD& b)
{	assert(a.ncols==b.nrows);
	MatrixD tmp(a.nrows,b.ncols);
	double* ptr;
	for(int i=0;i<a.nrows;i++)
		for(int j=0;j<b.ncols;j++)
		{	ptr=tmp.mat[i]+j;
			*ptr=0;
			for(int k=0;k<a.ncols;k++)
				(*ptr)+=a.mat[i][k]*b.mat[k][j];
		}
	return tmp;
}

ostream& operator << (ostream& out, const MatrixD& m)
{       for(int i=0;i<m.nrows;i++)
        {       out << "[ ";
                for(int j=0;j<m.ncols-1;j++) out << m.mat[i][j] << "  ";
                if(m.ncols>0) out << m.mat[i][m.ncols-1] << " ";
                out << "]" << endl;
        }
        return out;
}  

#endif