#include "Matrix.h"

using namespace std;

namespace numeric {

double ** Matrix::matrixmultiply( const Matrix & a, const Matrix & b)
{
	double ** result=0;
	if( a.dimCol != b.dimRow )
	{
		throw std::out_of_range("Dimension doesn't match");
	}
	else
	{
		result = _allocate(a.dimRow,b.dimCol);;
		for(unsigned i=0; i < a.dimRow; i++)
		{
			for(unsigned j=0; j < b.dimCol; j++)
			{
				for(unsigned k=0; k < a.dimCol; k++)
				{
            	result[i][j] += a(i,k)*b(k,j);
				}
			}
		}
	}
	return result;
}

void Matrix::_free(double **ptr,unsigned dim)
{
	if(ptr)
	{
		for(unsigned i=0;i<dim;i++)
			delete [] ptr[i];
		delete ptr;
	}
}

double ** Matrix::_allocate( unsigned rows,unsigned cols)
{
 	double ** ptr= new double*[rows];
 	for(unsigned i=0;i<rows;i++)
	{
		ptr[i] = new double[cols];
		for(unsigned j=0;j<cols;j++)
			ptr[i][j] = 0.0f;
	}
	return ptr;
}

double ** Matrix::_allocate( unsigned rows,unsigned cols,const Matrix & m)
{
 	double ** ptr= new double*[rows];
 	for(unsigned i=0;i<rows;i++)
	{
		ptr[i] = new double[cols];
		for(unsigned j=0;j<cols;j++)
			ptr[i][j] = m(i,j);
	}
	return ptr;
}

Matrix::Matrix()
:vector(0),
 dimRow(0),
 dimCol(0)
{}

Matrix::Matrix(unsigned rows,unsigned cols)
:dimRow(rows),
 dimCol(cols)
{
	vector = _allocate(dimRow,dimCol);
}

Matrix::Matrix(unsigned n)
:dimRow(n),
 dimCol(n)
{
	vector = _allocate(dimRow,dimCol);
}

Matrix::Matrix(const Matrix & m)
:vector(0),
 dimRow(0),
 dimCol(0)
{
		*this=m;
}

Matrix & Matrix::operator=(const Matrix & m)
{
	if( this != &m )
	{
		_free(vector,dimRow);
		vector=_allocate(m.dimRow,m.dimCol,m);
        dimRow = m.dimRow;
		dimCol = m.dimCol;
	}
	return *this;
}

double Matrix::operator()(unsigned row, unsigned col) const
{
	if( row < dimRow && col < dimCol )
		return vector[row][col];
	else
	{
		throw std::out_of_range("Dimension doesn't match in operator()");
		return 0.0;
	}
}

double & Matrix::operator()(unsigned row, unsigned col)
{
	if( row < dimRow && col < dimCol )
		return vector[row][col];
	else
	{
		throw std::out_of_range("Dimension doesn't match in operator()");
		return **vector;
	}
}

Matrix & Matrix::rightProduct( const Matrix & other)
{
	if( other.dimRow != dimCol )
	{
		throw std::out_of_range("Dimension doesn't match in rightProduct");
	}
	else
	{
		double **result = matrixmultiply(*this,other);
		_free(vector,dimRow);
		vector = result;
		dimCol = other.dimCol;
	}
	return (*this);
}

Matrix & Matrix::leftProduct(const Matrix & other)
{
	if( other.dimCol != dimRow )
	{
		throw std::out_of_range("Dimension doesn't match in leftProduct");
	}
	else
	{
		double **result = matrixmultiply(other,*this);
		_free(vector,dimRow);
		vector = result;
		dimRow = other.dimRow;
	}
	return (*this);
}

Matrix Matrix::transpose()
{
	Matrix result(dimCol,dimRow);
	for(unsigned i =0; i < dimCol; i++ )
		for(unsigned j =0; j < dimRow; j++ )
			 result.vector[i][j] = vector[j][i];
	return result;
}


std::ostream & operator<<( std::ostream & flux, const Matrix & a)
{
	for(unsigned i =0; i < a.getDimRow(); i++ )
	{
		for(unsigned j =0; j < a.getDimCol(); j++ )
		if( -1e-12 >= a(i,j) || a(i,j) > 1e-12)
			flux << a(i,j) << " ";
		else
		  flux << "0 " ;
		flux << std::endl;
	}

	return flux;
}

}