//written by Andre Betz 2007
//http://www.andrebetz.de

public class myMatrix {
	private int m_XSize;
	private int m_YSize;
	private double[] m_Field;
	private double m_Eps = 0.001;
	
	public myMatrix(int XSize,int YSize){
		m_XSize = XSize;
		m_YSize = YSize;
		m_Field = new double[m_XSize*m_YSize];
		for(int i=0;i<m_XSize*m_YSize;i++){
			m_Field[i] = 0.0;
		}
	}
	public myMatrix(myMatrix m){
		m_XSize = m.GetXSize();
		m_YSize = m.GetYSize();
		m_Field = new double[m_XSize*m_YSize];
		CopyFields(m.m_Field);
	}
	public myMatrix(double[][] fld){
		m_YSize = fld.length;
		m_XSize = fld[0].length;
		m_Field = new double[m_XSize*m_YSize];
		for(int y=0;y<m_YSize;y++){
			for(int x=0;x<m_XSize;x++){
				SetVal(x,y,fld[y][x]);
			}
		}
	}
	public int GetXSize(){
		return m_XSize;
	}
	public int GetYSize(){
		return m_YSize;
	}
	public double GetVal(int x,int y){
		double ret = 0.0;
		if(x<m_XSize&&y<m_YSize)
		{
			ret = m_Field[y*m_XSize+x];
		}
		return ret;
	}
	public double SetVal(int x,int y,double val){
		double ret = 0.0;
		if(x<m_XSize&&y<m_YSize)
		{
			m_Field[y*m_XSize+x] = val;
		}
		else
			throw new RuntimeException("Wrong Position");
		return ret;
	}
	public myMatrix Add(myMatrix m){
		if(m.GetXSize()!=m_XSize || m.GetYSize()!=m_YSize)
			throw new RuntimeException("Different Matrix Size");
		
		myMatrix n = new myMatrix(m_XSize,m_YSize);
		for(int y=0;y<m_YSize;y++){
			for(int x=0;x<m_XSize;x++){
				n.SetVal(x,y,this.GetVal(x, y) + m.GetVal(x, y));
			}
		}
		return n;
	}
	public myMatrix Minus(myMatrix m){
		if(m.GetXSize()!=m_XSize || m.GetYSize()!=m_YSize)
			throw new RuntimeException("Different Matrix Size");
		
		myMatrix n = new myMatrix(m_XSize,m_YSize);
		for(int y=0;y<m_YSize;y++){
			for(int x=0;x<m_XSize;x++){
				n.SetVal(x,y,this.GetVal(x, y) - m.GetVal(x, y));
			}
		}
		return n;
	}
	public myMatrix Mul(myMatrix m){
		if(m_XSize!=m.GetYSize())
			throw new RuntimeException("Different Matrix Size");
		myMatrix n = new myMatrix(m.GetXSize(),m_YSize);
		for(int y=0;y<n.GetYSize();y++){
			for(int x=0;x<n.GetXSize();x++){
				double tmp = 0.0;
				for(int i=0;i<m_XSize;i++){
					tmp += this.GetVal(i, y) - m.GetVal(x, i);
				}
				n.SetVal(x,y,tmp);
			}
		}
		return n;
	}
	public void SwapRow(int a,int b){
		if(a<m_YSize&&b<m_YSize)
		{
			for(int i=0;i<m_XSize;i++){
				double tmp = GetVal(i,a);
				SetVal(i,a,GetVal(i,b));
				SetVal(i,b,tmp);
			}
		}else
			throw new RuntimeException("Wrong Positions");
	}
	public String getString() {
		String str = "";
		for(int y=0;y<m_YSize;y++){
			for(int x=0;x<m_XSize;x++){
				str += GetVal(x,y) + "\t";
			}
			str += "\n";
		}
		return str;
    }
	public myMatrix Solve_Gauss(myMatrix vec)
	{
		if(vec.GetXSize()!=1 || vec.GetYSize() != m_YSize || m_YSize  != m_XSize)
			throw new RuntimeException("Wrong Matrix Dimensions");

		myMatrix A = new myMatrix(this);
		myMatrix b = new myMatrix(vec);
        myMatrix r = new myMatrix(m_XSize,1);

		for(int r1=0;r1<m_YSize;r1++){
            int maxPiv = r1;
            double MaxVal = Math.abs(A.GetVal(r1,maxPiv));
            for(int r2=r1+1;r2<m_XSize;r2++){            	
                if(Math.abs(A.GetVal(r1,r2)) > Math.abs(A.GetVal(r1,maxPiv))){
                	maxPiv = r2;
                	MaxVal = Math.abs(A.GetVal(r1,r2));
                }
            }
            
            if(MaxVal<m_Eps)
            	throw new RuntimeException("singular Matrix");
            
            A.SwapRow(r1, maxPiv);
            b.SwapRow(r1, maxPiv);

            for(int r2=r1+1;r2<m_XSize;r2++){ 
                double val = A.GetVal(r1,r2) / A.GetVal(r1,r1);
                for(int i=r1+1;i<m_XSize;i++) {
                	double t = A.GetVal(i,r2) - A.GetVal(i,r1) * val;
                	A.SetVal(i,r2,t);
                }
                A.SetVal(r1,r2,0.0);
				val = b.GetVal(0,r2) - b.GetVal(0,r1) * A.GetVal(r1,r2) / A.GetVal(r1,r1);
            	b.SetVal(0,r2,val);
            }
        }

        for(int y=m_YSize-1;y>=0;y--){
        	double t = 0.0;
        	for(int x=y+1;x<m_XSize;x++){
                t += A.GetVal(x,y) * r.GetVal(x,0);
        	}
        	double val = (b.GetVal(0,y) - t) / A.GetVal(y,y);
        	r.SetVal(y,0,val);
        }
        return r;
	}
	private void CopyFields(double[] fld){
		for(int i=0;i<m_XSize*m_YSize;i++){
			m_Field[i] = fld[i];
		}
	}
	
	public static void main(String[] args) {
		double[][] d = { 	{ -1, 3, 0,-10 }, 
							{ 0, 0, -3, 0 }, 
							{ 1, -3, 0,11 },
							{ 1, 0, 0, 1} };
		double [][] e = {{5},{46.5},{-6},{-20}};
		
		myMatrix A = new myMatrix(d);
		myMatrix b = new myMatrix(e);
		
		myMatrix c = A.Solve_Gauss(b);
		
		System.out.println(A.getString());
		System.out.println(b.getString());
		System.out.println(c.getString());
        System.out.println();
	}
}