/*
 * Created on 08.04.2005
 *
 * TODO To change the template for this generated file go to
 * Window - Preferences - Java - Code Style - Code Templates
 */

/**
 * @author GerberDom
 *
 * TODO To change the template for this generated type comment go to
 * Window - Preferences - Java - Code Style - Code Templates
 */
/* LUZerlegung_Variante01.java                    by Christof Senn, March 2003 */
/* --> Zerlegung einer nxn-Matrix in eine untere und eine obere Dreiecksmatrix */
/* --> in Variante 01 werden die Elt l_ii=1 gesetzt, für i={1, 2, 3, .., n}    */



public class LU_Zerlegung_Variante01 {


	/* Matrixgroesse: */
	private static int n = 4;
	private static final int vorKomma = 3;
	private static final int nachKomma= 5;

	public static void main(String[] args) {	
		double A[][] = 
		{	{2.1756,  4.0231,   -2.1732,  5.1967},
	 		{-4.0231, 6.0,      0.0,      1.1973},
	 		{-1.0,    -5.2107,  1.1111,   0.0},
			{6.0235,  7.0,      0.0,      -4.1561}	};
	
		double L[][] = new double[n][n];
		double U[][] = new double[n][n];
		double Kontrolle[][] = new double[n][n];
		for(int i = 0; i < n; i++){
			for(int j = 0; j < n; j++){
				L[i][j] = 0.0;
	   			U[i][j] = 0.0;
	   			Kontrolle[i][j] = 0.0;
			}
		}
			
		/* Matrixfaktorisierung: Zerlegung von A in L und U */
		U[0][0] = A[0][0];
		L[0][0] = 1.0;
		
		for(int j = 1; j < n; j++){
			U[0][j] = A[0][j];
			L[j][0] = A[j][0] / U[0][0];
			L[j][j] = 1.0;
		}
		
		for(int i = 1; i < n-1; i++){
			U[i][i] = A[i][i] - sum(L,U,i,i,i);
			for(int j= i+1; j < n; j++){
					U[i][j] = A[i][j] - sum(L,U,i,i,j);
				L[j][i] = 1.0 / U[i][i] * (A[j][i] - sum(L,U,i,j,i));
			}
		}
		U[n-1][n-1] = A[n-1][n-1] - sum(L,U,n-1,n-1,n-1);
		
		/* Kontrolle: */
		for(int i = 0; i < n; i++){
			for(int j = 0; j < n; j++){
				for(int k = 0; k < n; k++){
					Kontrolle[i][j] += L[i][k] * U[k][j];
				}
			}
		}
		
		/* Ausgabe der Lösungsmatrizen auf die Konsole: */
		InOut.print("\n"+n+" x "+n+" - Matrix: A\n\n");
		for(int i = 0; i < n; i++){
			for(int j = 0; j < n; j++){
				InOut.print(A[i][j], vorKomma, nachKomma);
				InOut.print("\t");
			}
			InOut.print("\n");
		}
		InOut.print("\n\n");
		
		InOut.print("Kontroll-Matrix: A = LU\n\n");
		for(int i = 0; i < n; i++){
			for(int j = 0; j < n; j++){
				InOut.print(Kontrolle[i][j], vorKomma, nachKomma);
				InOut.print("\t");
			}
			InOut.print("\n");
		}
		InOut.print("\n\n");
		
		InOut.print("\nL (Matrix)\n\n");
		for(int i = 0; i < n; i++){
			for(int j = 0; j < n; j++){
				InOut.print(L[i][j], vorKomma, nachKomma);
				InOut.print("\t");
			}
			InOut.print("\n");
		}
		InOut.print("\n\n");
		
		InOut.print("\nU (Matrix)\n\n");
		for(int i = 0; i < n; i++){
			for(int j = 0; j < n; j++){
				InOut.print(U[i][j], vorKomma, nachKomma);
				InOut.print("\t");
			}
			InOut.print("\n");
		}
		InOut.print("\n\n");
	}
	
	
	public static double sum(double L[][], double U[][], int bis, int i, int j) {
		double sum = 0;
		for(int k = 0; k < bis; k++){
			sum += L[i][k] * U[k][j];			
		}
		return sum;
	}
}