import java.util.Random;

public class QuadMatrix extends Matrix{
    
    
    private Matrix LR;
    private boolean isModified_LR;
    private boolean isModified_EW;
	private double klEw;
	private Vektor klEv;

	//////////////////////////////
	//		Konstruktoren...	//
	//////////////////////////////
    public QuadMatrix(int n){
		super(n,n);
		isModified_LR = true;
		isModified_EW = true;

    }
	public QuadMatrix(QuadMatrix m){
		super((Matrix) m);
		isModified_LR=true;
		isModified_EW = true;
	}


	////////////////////////////////////////////////////////////////////
	//		Methoden fuer die LR-Zerlegung dieser Matrix			  //
	////////////////////////////////////////////////////////////////////
	/*
	*berechnet x aus der Gleichung
	*Ax=b
	*mit Hilfe der LR-Zerlegung
	*/
	public Vektor Ax_b(Vektor b){
		if(isModified_LR)
			zerlegungLR();
		Vektor tmp = computeLy(b);
		return computeRx(tmp);

	}

	/*
	*Selektor, der die untere Dreiecksmatrix L der LR-Zerlegung von A zurueck gibt
	*so dass gilt
	* A=LR;
	*/
    public Matrix getL(){
		if(isModified_LR)
			zerlegungLR();
		QuadMatrix l = new QuadMatrix(this.matrix.length);
		for(int i=0; i<this.matrix.length; i++)
		    for(int j=0; j<this.matrix[i].length; j++){
				if(i==j)
				    l.matrix[i][j] = 1.0;
				if(i>j)
				    l.matrix[i][j] = LR.matrix[i][j]; 
		    }
		return l;
    }

	/*
	*gibt die obere Dreiecksmatrix R der LR-Zerlegung von A zurueck
	*so dass gilt
	* A=LR;
	*/
    public Matrix getR(){
		if(isModified_LR)
			zerlegungLR();
		QuadMatrix r = new QuadMatrix(this.matrix.length);
		for(int i=0; i<this.matrix.length; i++)
		    for(int j=0; j<this.matrix[i].length; j++){
				if(i<=j)
				    r.matrix[i][j] = LR.matrix[i][j]; 
			}
		return r;
    }



	/*
	*Die Implementierung des in der Vorlesung vorgestellten 
	*Algorithmus zur LR-Zerlegung.
	*ACHTUNG: ohne Spaltenpivotsuche!
	*/
    private void zerlegungLR(){
		if(isModified_LR){   
			//Methode arbeitet auf einer Kopie von A, damit A nicht zerstoert wird
		    LR = new QuadMatrix(this);
		
		    //////////////////////////////////////////////////////////
		    // hier beginnt die LR-Zerlegung ohne Spaltenpivotsuche!//
		    //////////////////////////////////////////////////////////
			 for(int k=0; k<(LR.matrix.length-1); k++)
				for(int i=k+1; i<LR.matrix.length; i++){
				    double l = LR.matrix[i][k]/LR.matrix[k][k];
				    LR.matrix[i][k] = l;
				    for(int j=k+1; j<LR.matrix.length; j++)
						LR.matrix[i][j]=LR.matrix[i][j]-l*LR.matrix[k][j];
				}
		    isModified_LR = false;   
		}
		return;
    }
    
    private Vektor computeLy(Vektor b){
		if(isModified_LR)
		   zerlegungLR();
		Vektor y = new Vektor(matrix.length);	
		for(int i=0;i<matrix.length;i++){
			double sum=0;
			for(int j=0;j<i;j++)//summe der bereits bekannten Produkte Xij*Rij
				sum = sum+y.getElement(j)*LR.matrix[i][j];
			y.setElement(b.getElement(i)-sum,i);
		}
		return y;
	
    }

 	
    private Vektor computeRx(Vektor y){
		if(isModified_LR)
		    zerlegungLR();
		Vektor x = new Vektor(matrix.length);	
		for(int i=0;i<matrix.length;i++){
			double sum=0;
			for(int j=0;j<i;j++)//summe der bereits bekannten Produkte Xij*Rij
				sum = sum+x.getElement(matrix.length-1-j)*LR.matrix[(matrix.length-1)-i][(matrix.length-1)-j];
			x.setElement((y.getElement((matrix.length-1)-i)-sum)/LR.matrix[(matrix.length-1)-i][(matrix.length-1)-i],(matrix.length-1)-i);
		}
		return x;
    }
    


	////////////////////////////////////////////////////
	//			Der Eigenwertalgorithmus			  //
	////////////////////////////////////////////////////

	/*
	*Selektor, um den kleinsten Eigenvektor zu erhalten
	*/
	public Vektor getKlEv(){
		if(isModified_EW)
			ew();

		return klEv;
	}


	/*
	*Selektor, um den kleinsten Eigenwert zu erhalten
	*/
	public double getKlEw(){
		if(isModified_EW)
			ew();
		return klEw;
	}
	

	/*
	*Hilfsmethde, um den sog. Raylay-Quotienten zu berechnen
	*/
	private double raylayQuot(Vektor x){
		double res = (Vektor.euklSkalarprodukt(x,Ax(x))/Vektor.euklSkalarprodukt(x,x));
		return res;
	}

	/*
	*Implementiert den in der Aufgabe beschrienbenen Algorithmus
	*zur bestimmung des kleinsten Eigenwerts und zugehoerigen Eigenvektors
	*/
	private void ew(){
		if(!isModified_EW)
			return;
		//pseudozufaelligen Vektor erzeugen
		klEv = new Vektor(matrix.length);
		Random r = new Random();
		do{
			for(int i=0;i<klEv.getLength();i++)
				klEv.setElement(r.nextDouble(),i);
		}while(Vektor.euklNorm(klEv)==0);
		
		if(!isModified_LR)
			zerlegungLR();

		//Initialisierung des Ergebnisses in Abhaengigkeit des Starvektors
		klEw = raylayQuot(klEv);
		while(Vektor.distance(Ax(klEv),klEv.mult(klEw))>0.000001){
			klEv=Ax_b(klEv);
			klEv=klEv.mult(1.0/Vektor.euklNorm(klEv));
			klEw = raylayQuot(klEv);
		}
		
	}


	////////////////////////////////////////////////////////////////////////
	//			allgemein benoetigte Operationen auf quadr. Matrizen	  //
	////////////////////////////////////////////////////////////////////////
	/*
	*Dies schreibt das Element x
	*in die (n,m)-te Komponente
	*/
	public void setElement(double x, int n, int m){
		super.setElement(x,n,m);
		isModified_LR=true;
		isModified_EW = true;
		return;
	}

	/*
	*gibt die Determinante zurueck
	*/
	public double detA(){
		if(isModified_LR)
			zerlegungLR();
		double det = 1;
		for(int i=0;i<matrix.length; i++)
			det=det*matrix[i][i];
		return det;
	}


}