/*
author:
	seth from http://www.bierdatenbank.de
tab size:
	2 (otherwise hardly readable!)
*/
#include "matrix.h"
using std::vector;

// zugriff per klammern ()
template<typename T> inline T& matrix<T>::operator() (unsigned row, unsigned col){
	if(row >= nrows_ || col >= ncols_) throw BoundsViolation();
	return data_[row*ncols_ + col];
}
template<typename T> inline const T& matrix<T>::operator() (unsigned row, unsigned col) const{
	if(row >= nrows_ || col >= ncols_) throw BoundsViolation();
	return data_[row*ncols_ + col];
}
// ctor
template<typename T> inline matrix<T>::matrix(unsigned nrows, unsigned ncols)
		: nrows_ (nrows), ncols_ (ncols){
	if(!nrows || !ncols) throw BadSize();
	data_ = new T[nrows * ncols];
	seth_std::mem_usage::used_mem(nrows*ncols * sizeof(T));
}
// copy-ctor
template<typename T> inline matrix<T>::matrix(const matrix<T>& m)
		: nrows_ (m.nrows_), ncols_ (m.ncols_){
	data_ = new T[nrows_ * ncols_];
	seth_std::mem_usage::used_mem(nrows_*ncols_ * sizeof(T));
	for(unsigned j,i=0;i<nrows_;++i)
		for(j=0;j<ncols_;++j)
			(*this)(i,j)=m(i,j);
}
// copy-ctor: import aus array
template<typename T> inline matrix<T>::matrix(T** array, unsigned nrows, unsigned ncols)
		: nrows_ (nrows), ncols_ (ncols){
	data_ = new T[nrows * ncols];
	seth_std::mem_usage::used_mem(nrows*ncols * sizeof(T));
	for(unsigned j,i=0;i<nrows;++i)
		for(j=0;j<ncols;++j)
			(*this)(i,j)=array[i][j];
}
// dtor
template<typename T> inline matrix<T>::~matrix(){
	delete[] data_;
	seth_std::mem_usage::released_mem(nrows_*ncols_ * sizeof(T));
}
// m = m (zuweisung)
template<typename T> inline matrix<T>& matrix<T>::operator=(const matrix<T>& m){
	if(nrows_==m.size(0) && ncols_==m.size(1)){
		for(unsigned j,i=0;i<nrows_;++i)
			for(j=0;j<ncols_;++j)
				(*this)(i,j)=m(i,j);
	}else std::cout <<"error: dimensionen bei addition muessen gleich sein."<<std::endl;
	return *this;
}

// -m : unitaerer matrix-operator; negiere komponentenweise
template<typename T> inline matrix<T> matrix<T>::operator-(){
	matrix<T> H(nrows_, ncols_);
	for(unsigned j,i=0;i<nrows_;++i)
		for(j=0;j<ncols_;++j)
			H(i,j)=-(*this)(i,j);
	return H;
}
// m*m = m
template<typename T> inline matrix<T> matrix<T>::operator*(const matrix<T>& m) const{
	matrix<T> H(nrows_, m.size_cols());
	if(ncols_==m.size_rows()){
		for(unsigned k,j,i=0;i<nrows_;++i)
			for(j=0;j<m.size_cols();++j){
				H(i,j)=0;
				for(k=0;k<ncols_;++k) H(i,j)+=(*this)(i,k)*m(k,j);
			}
	}else{
		std::cout <<"error: bei A*B muss gelten: spalten(A)==zeilen(B)"<<std::endl;
		H.zeros();
	}
	return H;
}
// m*v = v
template<typename T> inline vector<T> matrix<T>::operator*(const vector<T>& v) const{
	vector<T> v_result(nrows_);
	if(ncols_==v.size()){
		for(unsigned k,i=0;i<nrows_;++i){
			v_result.at(i)=0;
			for(k=0;k<ncols_;++k) v_result.at(i)+=(*this)(i,k)*v.at(k);
		}
	}else std::cout <<"error: bei M*v muss gelten: spalten(M)==zeilen(v)"<<std::endl;
	return v_result;
}
// m*s = m (komponentenweise)
template<typename T> inline matrix<T> matrix<T>::operator*(T skalar) const{
	matrix<T> H(*this);
	for(unsigned j,i=0;i<nrows_;++i)
		for(j=0;j<ncols_;++j)
			H(i,j)*=skalar;
	return H;
}
// m*=m = m
template<typename T> inline matrix<T>& matrix<T>::operator*=(const matrix<T>& m){
	*this=(*this)*m;
	return *this;
}
// m*=s = m
template<typename T> inline matrix<T>& matrix<T>::operator*=(T skalar){
	for(unsigned j,i=0;i<nrows_;++i)
		for(j=0;j<ncols_;++j)
			(*this)(i,j)*=skalar;
	return *this;
}
// m+m = m
template<typename T> inline matrix<T> matrix<T>::operator+(const matrix<T>& m) const{
	matrix<T> H(nrows_, ncols_);
	if(nrows_==m.size(0) && ncols_==m.size(1)){
		for(unsigned j,i=0;i<nrows_;++i)
			for(j=0;j<ncols_;++j)
				H(i,j)=(*this)(i,j)+m(i,j);
	}else{
		std::cout <<"error: dimensionen bei addition muessen gleich sein."<<std::endl;
		H.zeros();
	}
	return H;
}
// m+=m = m
template<typename T> inline matrix<T>& matrix<T>::operator+=(const matrix<T>& m){
	if(nrows_==m.size(0) && ncols_==m.size(1)){
		for(unsigned j,i=0;i<nrows_;++i)
			for(j=0;j<ncols_;++j)
				(*this)(i,j)+=m(i,j);
	}else std::cout <<"error: dimensionen bei addition muessen gleich sein."<<std::endl;
	return *this;
}
// m+m = m
template<typename T> inline matrix<T> matrix<T>::operator-(const matrix<T>& m) const{
	matrix<T> H(nrows_, ncols_);
	if(nrows_==m.size(0) && ncols_==m.size(1)){
		for(unsigned j,i=0;i<nrows_;++i)
			for(j=0;j<ncols_;++j)
				H(i,j)=(*this)(i,j)-m(i,j);
	}else{
		std::cout <<"error: dimensionen bei addition muessen gleich sein."<<std::endl;
		H.zeros();
	}
	return H;
}
// m+=m = m
template<typename T> inline matrix<T>& matrix<T>::operator-=(const matrix<T>& m){
	if(nrows_==m.size(0) && ncols_==m.size(1)){
		for(unsigned j,i=0;i<nrows_;++i)
			for(j=0;j<ncols_;++j)
				(*this)(i,j)-=m(i,j);
	}else std::cout <<"error: dimensionen bei addition muessen gleich sein."<<std::endl;
	return *this;
}
template<typename T> inline bool matrix<T>::operator==(const matrix<T>& m) const{
	bool result=(nrows_==m.size(0) && ncols_==m.size(1));
	if(result)
		for(unsigned j,i=0;result && i<nrows_;++i)
			for(j=0;result && j<ncols_;++j)
				result=((*this)(i,j)==m(i,j));
	return result;
}
template<typename T> inline bool matrix<T>::operator!=(const matrix<T>& m) const{
	return !((*this)==m);
}

// ausgabe
std::ostream& operator<<(std::ostream& strm, const matrix<int>& m){
	for(unsigned j,i=0;i<m.size_rows();++i){
		strm <<"[";
		for(j=0;j<m.size_cols()-1;++j)
			strm <<m(i,j)<<" ";
		strm <<m(i,j)<<"]"<<std::endl;
	}
	return strm;
}
std::ostream& operator<<(std::ostream& strm, const matrix<double>& m){
	for(unsigned j,i=0;i<m.size_rows();++i){
		strm <<"[";
		for(j=0;j<m.size_cols()-1;++j)
			strm <<m(i,j)<<" ";
		strm <<m(i,j)<<"]"<<std::endl;
	}
	return strm;
}
// ueberpruefe, ob dimensionen eingehalten
template<typename T> inline bool matrix<T>::isIndex(unsigned i, unsigned j) const{
	return (i<nrows_ && j<ncols_);
}
// anzahl zeilen(0) oder spalten(1
template<typename T> inline unsigned matrix<T>::size(bool cols) const{
	return (cols)?size_cols():size_rows();
}
// anzahl zeilen
template<typename T> inline unsigned matrix<T>::size_rows() const{
	return nrows_;
}
// anzahl spalten
template<typename T> inline unsigned matrix<T>::size_cols() const{
	return ncols_;
}
// aendert matrix zu identitaet
template<typename T> inline void matrix<T>::identity(){
	for(unsigned j,i=0;i<nrows_;++i)
		for(j=0;j<ncols_;++j)
			(*this)(i,j)=(i==j);
}
// invertiert nxn-matrix mit gauss-jordan-algorithmus mit spalten-pivotisierung
template<typename T> inline matrix<T> matrix<T>::inverse() const{
	matrix<T> H(*this);
	if(nrows_==ncols_){
		unsigned col,row,j,el=0;
		unsigned *v;
		v=new unsigned[nrows_];
		T b,max_el;
		for(col=0;col<nrows_;++col){
			v[col]=col;
			max_el=0;
			for(row=col;row<nrows_;++row)	// suche max_el:=betragsgroesstes elem. in spalte col, ab col-tem element
				if(fabs(H(row,col))>fabs(max_el)){
					max_el=H(row,col);
					el=row;										// max_el = el-te komponente in col-ter spalte
				}
			if(max_el==0){								// (nullspalte =>) matrix singulaer
				std::cout <<"error: matrix nicht invbar"<<std::endl;
				return (*this);
			}
			if(el!=col){									// wenn max_el nicht ganz oben (in col-ter colxcol-hauptuntermatrix) steht vertausche zeilen, so dass max_el ganz oben steht also =a[col][col]
				for(j=0;j<nrows_;++j) std::swap(H(col,j), H(j,el));
				v[col]=el;									// speichere nummer der gewechselten zeile
			}
			H(col,col)=1;									// setze max_el auf 1
			for(j=0;j<nrows_;++j) H(col,j)/=max_el;// normalisiere zeile col (mit 1/max_el)
			for(row=0;row<nrows_;++row)
				if(row!=col){								// in allen anderen zeilen
					b=H(row,col);
					H(row,col)=0;							// anulliere spalte col
					for(j=0;j<nrows_;++j) H(row,j)-=b*H(col,j);
				}
		}
		for(col=nrows_-1;col>=0;--col){
			if(v[col]!=col)
				for(j=0;j<nrows_;++j) std::swap(H(j,col), H(j,v[col]));
			if(col==0) break;
		}
		delete[] v;
	}else std::cout <<"error: matrix nicht quadratisch => nicht invbar"<<std::endl;
	return H;
}
// loese Ax=b
template<typename T> inline void matrix<T>::solve(const vector<T>& b, vector<T>& x) const{
	x=(*this).inverse()*b;
}
// aendert matrix zu nullmatrix
template<typename T> inline void matrix<T>::zeros(){
	for(unsigned j,i=0;i<nrows_;++i)
		for(j=0;j<ncols_;++j)
			(*this)(i,j)=0;
}
// transponiert matrix
template<typename T> inline matrix<T> matrix<T>::transpose(){
	matrix H(ncols_, nrows_); // ncols und nrows vertauscht
	for(unsigned j,i=0;i<nrows_;++i)
		for(j=0;j<ncols_;++j)
			H(j,i)=(*this)(i,j);
	return H;
}

// damit die linker-errors verschwinden
template matrix<seth_std::punkt2d<int> >;
template matrix<seth_std::punkt2d<double> >;
template matrix<double>;
template matrix<int>;