package dmrgvisu2;

import org.netlib.blas.DGEMM;
import org.netlib.lapack.DGESVD;
import org.netlib.lapack.DSYEV;
import org.netlib.util.*;

public class Matrix
{
	private static final long serialVersionUID = 1L;
	int [] Size;	
	double[][] Content;
/************************************************************
*	Usage:													*
*	#========== x =====================================>	*
*	|[0,0]	[0,1]	[0,2]	...	...	[0,Size[1]]				*
*	|[1,0]								.					*
*	|   .								.					*
*	y  .								.					*
*	|   .								.					*
*	|[Size[0],0]	...	...	...	[Size[0],Size[1]]			*
*	v														*
*															*
************************************************************/
	int[] Merges;
	boolean[] MergesUnder;

	public int[] GetSize()
	{
		return this.Size;
	}
	public double[][] GetContent()
	{
		return this.Content;
	}

	public Matrix()
	{
		this.SetSize(1,1);
		this.Merges = new int[0];
		this.MergesUnder = new boolean[0];
	}
	public Matrix(int[] Size)
	{
		this.SetSize(Size);
		this.Merges = new int[0];
		this.MergesUnder = new boolean[0];
	}
	public Matrix(int x, int y)
	{
		this.SetSize(x,y);
		this.Merges = new int[0];
		this.MergesUnder = new boolean[0];
	}
	public Matrix(Matrix NewMatrix)
	{
		this.Size = NewMatrix.GetSize();
		this.Content = NewMatrix.GetContent();
		this.Merges = NewMatrix.Merges;
		this.MergesUnder = NewMatrix.MergesUnder;
	}
	private void CreateContent()
	{
		this.Content = new double[this.Size[0]][this.Size[1]];
		for (int y = 0; y < this.Size[0]; y++)
			for (int x = 0; x < this.Size[1]; x++)
				this.Content[y][x] = 0;
	}
	public void SetContent(int x, int y, double Value)
	{
		this.Content[y][x] = Value;
	}
	public double GetContent(int x, int y)
	{
		return this.Content[y][x];
	}
	public void SetContent(double[][] Content)
	{
		this.SetSize(Content[0].length, Content.length);
		this.Content = Content;
	}
	public void SetSize(int x, int y)
	{
		int[] newSize = new int[2];
		newSize[0] = y;
		newSize[1] = x;
		this.SetSize(newSize);
	}
	public void SetSize(int[] Size)
	{
		this.Size = Size;
		this.CreateContent();
	}

	public String toString()
	{
		String newString = "This Matrix has the size: " + this.Size[1] + " by " + this.Size[0];
		if (this.Merges.length !=0)
		{
			newString = newString + " and is merged on ";
			for (int i=0; i<this.Merges.length;i++)
				newString = newString + this.Merges[i] + "\n";
		}
		else
			newString = newString + "\n";
		for (int y = 0; y < this.Size[0]; y++)
		{
			for (int x = 0; x < this.Size[1]; x++)
			{
				newString = newString + " " + this.Content[y][x];
			}
			newString = newString + "\n";
		}
		return newString;

	}
	
	public void printreadable()
	{
		System.out.format("This Matrix has the size: %d by %d", this.Size[1], this.Size[0]);
		if (this.Merges.length !=0)
		{
			System.out.format(" and is merged on ");
			for (int i=0; i<this.Merges.length;i++)
				System.out.format("%d %n", this.Merges[i]);
		}
		else
			System.out.format("%n");
		for (int y = 0; y < this.Size[0]; y++)
		{
			for (int x = 0; x < this.Size[1]; x++)
			{
				System.out.format(" %12.3f", this.Content[y][x]);
			}
			System.out.format("%n");
		}
		return;
	}
	
	public Matrix mult(Matrix b)
	{
		if (b.GetSize()[0]==1 && b.GetSize()[1]==1)
			return this.mult(b.GetContent(0,0));
		else if (this.GetSize()[0]==1 && this.GetSize()[1]==1)
			return b.mult(this.GetContent(0,0));
		else
			return this.mult(b, "n");
	}
	public Matrix mult(Matrix b, String transb)
	{
		String transa = "n";
		int m = (transa == "n"?this.GetSize()[0]:this.GetSize()[1]);
		int n = (transb == "n"?b.GetSize()[1]:b.GetSize()[0]);
		int k = this.GetSize()[1];
		double alpha = 1.0;
		double beta = 0.0;

		double [][] c = new double[m][n];
		DGEMM.DGEMM(transa, transb, m, n, k, alpha, this.GetContent(), b.GetContent(), beta, c);

		Matrix newMatrix = new Matrix(m,n);
		newMatrix.SetContent(c);
		return newMatrix;
	}
	public Matrix trans()
	{
		Matrix l = new Matrix(this.GetSize()[1], this.GetSize()[1]);
		for (int i = 0; i < this.GetSize()[1]; i++)
			l.SetContent(i,i,1);
		l = l.mult(this, "t");
		l.Merges = this.Merges;
		l.MergesUnder = this.MergesUnder;
		for (int i=0; i<this.MergesUnder.length; i++)
			l.MergesUnder[i] = !this.MergesUnder[i];
		return (l);
	}
	public Matrix add(Matrix b)
	{
		if ((this.GetSize()[0] != b.GetSize()[0]) || (this.GetSize()[1] != b.GetSize()[1]))
		{	
			System.out.println("ERROR: Dimension passt nicht");
		}
		Matrix newMatrix = new Matrix(this.GetSize());
		for (int x = 0; x < this.GetSize()[1]; x++)
			for (int y = 0; y < this.GetSize()[0]; y++)
				newMatrix.SetContent(x, y, this.GetContent(x,y) + b.GetContent(x,y));
		return newMatrix;
	}
	public Matrix mult(double d)
	{
		Matrix newMatrix = new Matrix(this.GetSize());
		for (int x = 0; x < this.GetSize()[1]; x++)
			for (int y = 0; y < this.GetSize()[0]; y++)
				newMatrix.SetContent(x, y, this.GetContent(x,y) * d);
		return newMatrix;
	}
	public Matrix neg()
	{
		Matrix newMatrix = new Matrix(this.GetSize());
		for (int x = 0; x < this.GetSize()[1]; x++)
			for (int y = 0; y < this.GetSize()[0]; y++)
				newMatrix.SetContent(x, y, -this.GetContent(x,y));
		return newMatrix;
	}
	public Matrix sub(Matrix b)
	{
		return this.add(b.neg());
	}
	public Matrix svd(Matrix u, Matrix vt)
	{
//		System.out.println("HIER IST DIE SVD:");
//		System.out.println(this);
		String jobU="S", jobVt="S";
		double[] s = new double[(this.Size[1]<this.Size[0]?this.Size[1]:this.Size[0])];
		int lwork = 5*s.length+(this.Size[1]>this.Size[0]?this.Size[1]:this.Size[0]);
	    double []work = new double[lwork];
		intW info = new intW(0);
		
		if (this.Size[1] < this.Size[0])
		{
			u.SetSize(this.Size[1], this.Size[0]);
			vt.SetSize(this.Size[1], this.Size[1]);
		}
 		else
		{
			u.SetSize(this.Size[0], this.Size[0]);
			vt.SetSize(this.Size[1], this.Size[0]);
		}
		DGESVD.DGESVD(jobU, jobVt, this.Size[0], this.Size[1], this.GetContent(), s, u.GetContent(), vt.GetContent(), work, lwork, info);

		u.Merges = this.Merges;
		u.MergesUnder = this.MergesUnder;
		vt.Merges = this.Merges;
		vt.MergesUnder = this.MergesUnder;

		Matrix sigma=new Matrix(s.length, s.length);
		for (int i=0; i<s.length; i++)
			sigma.SetContent(i, i, s[i]);

		return sigma;
	}

	public int merge(Matrix b, boolean under)
	{
		if (under == true)
		{
			if (b.GetSize()[1] != this.GetSize()[1])
				return -1;
//			System.out.println("this.GetSize()[0]: " + this.GetSize()[0] + ", this.GetSize()[1]: " + this.GetSize()[1] + ", b.GetSize()[0]: " + b.GetSize()[0] + ", b.GetSize()[1]: " + b.GetSize()[1]);
			Matrix newMatrix = new Matrix(this.GetSize()[1], this.GetSize()[0]+b.GetSize()[0]);
			newMatrix.Merges = new int[this.Merges.length+1];
			newMatrix.Merges[this.Merges.length] = this.GetSize()[0];
			newMatrix.MergesUnder = new boolean[this.MergesUnder.length+1];
			newMatrix.MergesUnder[this.MergesUnder.length] = true;
			for (int y = 0; y < this.GetSize()[0]; y++)
				for (int x = 0; x < this.GetSize()[1]; x++)
					newMatrix.SetContent(x,y,this.GetContent(x,y));
			for (int y = 0; y < b.GetSize()[0]; y++)
				for (int x = 0; x < b.GetSize()[1]; x++)
					newMatrix.SetContent(x,y+this.GetSize()[0], b.GetContent(x,y));
			this.Content = newMatrix.GetContent();
			this.Size = newMatrix.GetSize();
			this.Merges = newMatrix.Merges;
			this.MergesUnder = newMatrix.MergesUnder;
		}
		else
		{
			if (b.GetSize()[0] != this.GetSize()[0])
				return -1;
//			System.out.println("this.GetSize()[0]: " + this.GetSize()[0] + ", this.GetSize()[1]: " + this.GetSize()[1] + ", b.GetSize()[0]: " + b.GetSize()[0] + ", b.GetSize()[1]: " + b.GetSize()[1]);
			Matrix newMatrix = new Matrix(this.GetSize()[1]+b.GetSize()[1], this.GetSize()[0]);
			newMatrix.Merges = new int[this.Merges.length+1];
			newMatrix.Merges[this.Merges.length] = this.GetSize()[1];
			newMatrix.MergesUnder = new boolean[this.MergesUnder.length+1];
			newMatrix.MergesUnder[this.MergesUnder.length] = false;
			for (int y = 0; y < this.GetSize()[0]; y++)
				for (int x = 0; x < this.GetSize()[1]; x++)
					newMatrix.SetContent(x,y, this.GetContent(x,y));
			for (int y = 0; y < b.GetSize()[0]; y++)
				for (int x = 0; x < b.GetSize()[1]; x++)
					newMatrix.SetContent(x+this.GetSize()[1],y, b.GetContent(x,y));
			this.Content = newMatrix.GetContent();
			this.Size = newMatrix.GetSize();
			this.Merges = newMatrix.Merges;
			this.MergesUnder = newMatrix.MergesUnder;
		}
		return 0;
	}
	
	public Matrix split()
	{
		if (this.Merges.length == 0)
			return this;
		if (this.MergesUnder[this.MergesUnder.length-1] == true)
		{
			Matrix newMatrix = new Matrix(this.GetSize()[1], this.Merges[this.Merges.length-1]);
			Matrix outMatrix = new Matrix(this.GetSize()[1], this.GetSize()[0]-this.Merges[this.Merges.length-1]);
			
			newMatrix.Merges = new int[this.Merges.length-1];
			for (int i=0;i<this.Merges.length-1;i++)
				newMatrix.Merges[i] = this.Merges[i];
			newMatrix.MergesUnder = new boolean[this.MergesUnder.length-1];
			for (int i=0;i<this.MergesUnder.length-1;i++)
				newMatrix.MergesUnder[i] = this.MergesUnder[i];			
				
			for (int y = 0; y < newMatrix.GetSize()[0]; y++)
				for (int x = 0; x < newMatrix.GetSize()[1]; x++)
					newMatrix.SetContent(x,y, this.GetContent(x,y));
			for (int y = 0; y < outMatrix.GetSize()[0]; y++)
				for (int x = 0; x < outMatrix.GetSize()[1]; x++)
					outMatrix.SetContent(x,y, this.GetContent(x,y+newMatrix.GetSize()[0]));
			
			this.Content = newMatrix.GetContent();
			this.Size = newMatrix.GetSize();
			this.Merges = newMatrix.Merges;
			this.MergesUnder = newMatrix.MergesUnder;
//			System.out.println("HIER IST SPLITTING!");
//			System.out.println(outMatrix);
			return outMatrix;
		}
		else
		{
			Matrix newMatrix = new Matrix(this.Merges[this.Merges.length-1], this.GetSize()[0]);
			Matrix outMatrix = new Matrix(this.GetSize()[1]-this.Merges[this.Merges.length-1], this.GetSize()[0]);
			
			newMatrix.Merges = new int[this.Merges.length-1];
			for (int i=0;i<this.Merges.length-1;i++)
				newMatrix.Merges[i] = this.Merges[i];
			newMatrix.MergesUnder = new boolean[this.MergesUnder.length-1];
			for (int i=0;i<this.MergesUnder.length-1;i++)
				newMatrix.MergesUnder[i] = this.MergesUnder[i];			
				
			for (int y = 0; y < newMatrix.GetSize()[0]; y++)
				for (int x = 0; x < newMatrix.GetSize()[1]; x++)
					newMatrix.SetContent(x,y, this.GetContent(x,y));
//			System.out.println("this.GetSize()[0]: " + this.GetSize()[0] + ", this.GetSize()[1]: " + this.GetSize()[1] + ", outMatrix.GetSize()[0]: " + outMatrix.GetSize()[0] + ", outMatrix.GetSize()[1]: " + outMatrix.GetSize()[1] + ", newMatrix.GetSize()[0]: " + newMatrix.GetSize()[0] + ", newMatrix.GetSize()[1]: " + newMatrix.GetSize()[1]);
			for (int y = 0; y < outMatrix.GetSize()[0]; y++)
				for (int x = 0; x < outMatrix.GetSize()[1]; x++)
					outMatrix.SetContent(x,y, this.GetContent(x+newMatrix.GetSize()[1],y));
			
			this.Content = newMatrix.GetContent();
			this.Size = newMatrix.GetSize();
			this.Merges = newMatrix.Merges;
			this.MergesUnder = newMatrix.MergesUnder;
			return outMatrix;
		}
	}
	
	public double Trace()
	{
		double out=0.0;
		if (this.GetSize()[0]!=this.GetSize()[1])
			return 0.0;
		for (int i=0; i<this.GetSize()[0]; i++)
			out += this.GetContent(i,i);
		return out; 
	}
	
	public double[] Eigen()
	{
		if (this.GetSize()[0] != this.GetSize()[1])
			return new double[0];
		double []eva =new double[this.GetSize()[0]];
	    int lwork = this.GetSize()[0]*this.GetSize()[0]+2;
	    double []work = new double[lwork];
	    intW info = new intW(0);
		DSYEV.DSYEV("V", "U", this.GetSize()[0], this.GetContent(), eva, work, lwork, info);
		return eva;
	}
}

