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

/**
 * @author Gerry Broennimann, ib03
 */

/*  LUZerlegung_Cholesky.java    by Gerry Brönnimann, April 2005 */
/*  based on LUZerlegung_Variante02.java                   by Christof Senn, March 2003 */
/* --> Zerlegung einer nxn-Matrix in eine untere und eine obere Dreiecksmatrix */
/* --> in Variante 02 werden die Elt u_ii=1 gesetzt, für i={1, 2, 3, .., n}    */

public class LU_Zerlegung_Cholesky{
	/* Matrixgroesse: */
	private static int n = 10;
	private static final int vorKomma = 2;
	private static final int nachKomma= 3;

	public static void main(String[] args) {


		double A[][] = new double[10][10];

		// Matrix initialisieren:
		for (int i= 1; i<=10; i++)
		{
			for (int k=1; k<=10; k++)
			{
				if (i == k) A[i-1][k-1] = 1;
				else A[i-1][k-1] = (i+k-1)*0.05;
			}
		}



		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 */
			L[0][0] = A[0][0];
			U[0][0] = 1.0;

		/*
		for(int j = 1; j < n; j++){
			L[j][0] = A[j][0];
			U[0][j] = A[0][j] / L[0][0];
			U[j][j] = 1.0;
		}

		for(int i = 1; i < n-1; i++){
			L[i][i] = A[i][i] - sum(L,U,i,i,i);
			for(int j= i+1; j < n; j++){
				L[j][i] = A[j][i] - sum(L,U,i,j,i);
				U[i][j] = 1.0 / L[i][i] * (A[i][j] - sum(L,U,i,i,j));
			}
		}
		L[n-1][n-1] = A[n-1][n-1] - sum(L,U,n-1,n-1,n-1); */

		// Cholesky:
		L[0][0] = Math.sqrt(A[0][0]);

		for(int j = 1; j < n; j++){
			L[j][0] = A[j][0] / L[0][0];
		}

		for(int j = 1; j < n; j++){
					//L[j][0] = A[j][0] / L[0][0];
					//U[0][j] = A[0][j] / L[0][0];
					//U[j][j] = 1.0;
				}

				for(int i = 1; i < n; i++){
					L[i][i] = Math.sqrt( A[i][i] - sum(L,L,i,i,i) );
					for(int j= i+1; j < n; j++){
						//L[j][i] = A[j][i] - sum(L,U,i,j,i);
						L[i][j] = 1.0 / L[i][i] * (A[i][j] - sum(L,L,i,i,j));
					}
				}
		L[n-1][n-1] = Math.sqrt( A[n-1][n-1] - sum(L,L,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");
	}


	private 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;
	}
}