#include <iostream>
#include <math.h>
#include "matrix.h"
using namespace std;

Matrix::Matrix(int d):cross(0){
	this->dimension=d;
	this->elements=new double*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new double[d];
		for(j=0;j<d;j++)
			this->elements[i][j]=0;
	}	
}

Matrix::Matrix(int d,double **es):cross(0){
	this->dimension=d;
	this->elements=new double*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new double[d];
		for(j=0;j<d;j++)
			this->elements[i][j]=es[i][j];
	}
}

Matrix::Matrix(int d,double e,int c):cross(0){ //0 for every element, 1 just for cross
	this->dimension=d;
	this->elements=new double*[d];
	int i,j;
	if(!c)
		for(i=0;i<d;i++){
			this->elements[i]=new double[d];
			for(j=0;j<d;j++)
				this->elements[i][j]=e;	
		}
	else{
		this->cross=1;
		for(i=0;i<d;i++){
			this->elements[i]=new double[d];
			for(j=0;j<d;j++)
				if(i==j)this->elements[i][j]=e;	
				else this->elements[i][j]=0;
		}
	}
}

Matrix::Matrix(int d,double *&es,Fill_Way fw):cross(0){
	if(fw!=crs&&fw!=row&&fw!=col)d=int(sqrt(d));	
	this->dimension=d;
	this->elements=new double*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new double[d];	
		switch(fw){
			case crs:
				this->cross=1;
				for(j=0;j<d;j++){
					if(i==j)this->elements[i][j]=es[i];
					else this->elements[i][j]=0;	
				}
				break;
			case row:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[j];
				}
				break;
			case col:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i];
				}
				break;
			case all:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i*d+j];
				}
				break;
			default:
				cerr << "Warning: unknown fill way, use all as default!" << endl;
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i*d+j];
				}		
		}				
	}
}

Matrix::Matrix(int d,double *es,Fill_Way fw):cross(0){
	if(fw!=crs&&fw!=row&&fw!=col)d=int(sqrt(d));	
	this->dimension=d;
	this->elements=new double*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new double[d];	
		switch(fw){
			case crs:
				this->cross=1;
				for(j=0;j<d;j++){
					if(i==j)this->elements[i][j]=es[i];
					else this->elements[i][j]=0;	
				}
				break;
			case row:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[j];
				}
				break;
			case col:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i];
				}
				break;
			case all:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i*d+j];
				}
				break;
			default:
				cerr << "Warning: unknown fill way, use all as default!" << endl;
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i*d+j];
				}		
		}				
	}
}

Matrix::Matrix(const Matrix &m1,const Matrix &m2,const Matrix &m3,const Matrix &m4):cross(0){	
	if(m1.dimension-m2.dimension||m1.dimension-m3.dimension||m2.dimension-m3.dimension||
	   m1.dimension-m4.dimension||m2.dimension-m4.dimension||m3.dimension-m4.dimension){
		cerr <<	"Warning: dimension of 4 matrix are not equal, the this matrix!" << endl;
		exit(1);
	}
	int d=m1.dimension;
	this->dimension=2*d;
	this->elements=new double*[2*d];
	int i,j;
	for(i=0;i<2*d;i++){
		this->elements[i]=new double[2*d];		
		for(j=0;j<2*d;j++){
			if(i<d){
				if(j<d)this->elements[i][j]=m1.elements[i][j];
				else this->elements[i][j]=m2.elements[i][j-d];
			}
			else{
				if(j<d)this->elements[i][j]=m3.elements[i-d][j];
				else{ this->elements[i][j]=m4.elements[i-d][j-d];		
				}
			}
		}
	}//this->elements[0][0].print();
	//this->setCross();	
}

int Matrix::getDimension(void){
	return this->dimension;
}

double Matrix::getElement(int i,int j){
	if(i>=this->dimension||j>=this->dimension){
		cerr << "Warning: subscript out of boundary, return the first element!" << endl;
		return this->elements[0][0];
	}
	return this->elements[i][j];
}

void Matrix::setElement(double e,int i,int j){
	if(i>=this->dimension||j>=this->dimension){
		cerr << "Warning: subscript out of boundary, no elements set!" << endl;
		return;
	}	
	this->elements[i][j]=e;
	this->setCross();
}

int Matrix::operator==(Matrix &m){
	if(this->dimension-m.dimension)return 0;
	int i,j;
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			if(this->elements[i][j]-m.getElement(i,j)==0)return 0;
		}
	}
	return 1;
}

Matrix Matrix::operator/(double l){
	Matrix *temp=new Matrix(this->dimension);
	int i,j;
	for(i=0;i<this->dimension;i++)
		for(j=0;j<this->dimension;j++)
			temp->elements[i][j]=this->elements[i][j]/l;
	return *temp;
}

void Matrix::LU(Matrix *&L,Matrix *&U){
	int i,j,m;
	double s1,s2;
	L=new Matrix(this->dimension,1.0,1);
	U=new Matrix(this->dimension);
	for(i=0;i<this->dimension;i++){
		U->elements[0][i]=this->elements[0][i];
		if(i>0)L->elements[i][0]=this->elements[i][0]/U->elements[0][0];
	}
	for(i=1;i<this->dimension;i++){
		for(j=i;j<this->dimension;j++){
			s1=0;s2=0;
			for(m=0;m<i;m++){
				s1+=L->elements[i][m]*U->elements[m][j];
				if(j>i)s2+=L->elements[j][m]*U->elements[m][i];
			}			
			U->elements[i][j]=this->elements[i][j]-s1;
			if(j>i)L->elements[j][i]=(this->elements[j][i]-s2)/U->elements[i][i];
		}
	}
}

void Matrix::exRow(int m,int n){
	double *t;
	t=this->getRow(m);
	this->setRow(this->getRow(n),m);
	this->setRow(t,n);
	delete t;
}

void Matrix::setCross(void){
	int i,j;
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			if(i-j){
				if(!(this->elements[i][j]==0)){
					this->cross=0;return;
				};
			}
		}
	}
	this->cross=1;
}

int Matrix::getCross(void){
	return this->cross;
}

double *Matrix::getRow(int i){
	if(i>=this->dimension){
		cerr << "Warning: subscript out of boundary, return the first row!" << endl;
		return this->elements[0];
	}
	return this->elements[i];
}

double *Matrix::getCol(int j){
	double *tempCol;
	int i;
	tempCol=new double[this->dimension];
	if(j>=this->dimension){
		cerr << "Warning: subscript out of boundary, return the first row!" << endl;
		return this->getCol(0);
	}
	for(i=0;i<this->dimension;i++){
		tempCol[i]=this->elements[i][j];
	}
	return tempCol;
}

void Matrix::setRow(double *es,int i){
	int j;
	if(i>=this->dimension){
		cerr << "Warning: subscript out of boundary, no rows set!" << endl;
		return;
	}
	for(j=0;j<this->dimension;j++)
		this->setElement(es[j],i,j);
	this->setCross();
}

void Matrix::setCol(double *es,int j){
	int i;
	if(j>=this->dimension){
		cerr << "Warning: subscript out of boundary, no colomns set!" << endl;
		return;
	}	
	for(i=0;i<this->dimension;i++)
		this->setElement(es[i],i,j);
	this->setCross();
}

double *Matrix::getCrossE(void){
	double *tempCross;
	int i,j;
	tempCross=new double[this->dimension];
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			if(i==j){
				tempCross[i]=this->elements[i][j];
				break;
			}
		}
	}
	return tempCross;
}

Matrix Matrix::trans(void){
	int i,j;
	Matrix *r=new Matrix(this->dimension);
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			r->setElement(this->elements[i][j],j,i);
		}
	}
	return *r;
}

double *Matrix::operator*(double *s){
	int i,j;
	double *temp=new double[this->dimension];
	for(i=0;i<this->dimension;i++){
		temp[i]=0;
		for(j=0;j<this->dimension;j++){
			temp[i]=temp[i]+(this->elements[i][j]*s[j]);
		}
	}
	return temp;
}

Matrix Matrix::operator*(double &l){
	Matrix *r=new Matrix(this->dimension);
	int i,j;
	for(i=0;i<this->dimension;i++)
		for(j=0;j<this->dimension;j++)
			r->elements[i][j]=l*this->elements[i][j];
	return *r;
}

Matrix Matrix::operator*(const Matrix &m){
	if(this->dimension-m.dimension){
		cerr << "Warning: dimensions of two matrix not equal, return the left matrix!" << endl;
		return *this;
	}
	int i,j;
	Matrix *r=new Matrix(this->dimension); 
	double *es,*cs;
	es=new double[this->dimension];
	cs=new double[this->dimension];
	if(this->cross){					//左乘乘行 
		cs=this->getCrossE();
		for(i=0;i<this->dimension;i++){
			//es=m.getRow(i);
			for(j=0;j<this->dimension;j++){
				es[j]=cs[i]*m.elements[i][j];
			}
			r->setRow(es,i);
		}
	}
	else if(m.cross){					//右乘乘列 
		//cs=m.getCrossE();
		for(i=0;i<this->dimension;i++){
			es=this->getCol(i);
			for(j=0;j<this->dimension;j++){
				if(i==j)es[j]=m.elements[i][j]*es[j];
			}
			r->setCol(es,i);
		}
	}
	else{								//直接乘法 
		int s;
		for(s=0;s<this->dimension;s++){
			for(i=0;i<this->dimension;i++){
				for(j=0;j<this->dimension;j++){
					r->setElement(r->getElement(i,j)+this->elements[i][s]*m.elements[s][j],i,j);
				}		
			}			
		}
	}
	r->setCross();
	return *r;
}

Matrix Matrix::operator-(Matrix &m){
	Matrix *r=new Matrix(this->dimension);
	int i,j;
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			r->setElement(this->elements[i][j]-m.getElement(i,j),i,j);
		}
	}
	r->setCross();
	return *r;
}

void Matrix::setCrossE(double *es){
	int i,j;	
	for(i=0;i<this->dimension;i++)
		for(j=0;j<this->dimension;j++)
			if(i==j)this->setElement(es[i],i,j);
	this->setCross();
}

void Matrix::print(int d){
	int i,j;
	if(d>this->dimension)this->print();
	else
		for(i=0;i<d;i++){
			for(j=0;j<d;j++){
				cout << this->elements[i][j];
				cout << "\t";
			}
			cout << endl;
		}
}

void Matrix::print(void){
	this->print(this->dimension);
}

Matrix::~Matrix(void){
	for(int i=0;i<this->dimension;i++)
		delete [] this->elements[i];
	delete [] this->elements;		
}