clear;

A = rand(100,100);
n = size(A,1); % Groesse der Matrix
b = rand(n,1); % rechte Seite

A_copy = A; % zur spaeteren Ueberpruefung sichern wir die Matrix A

for k = 1:n-1
    for i = k+1:n
        A(i,k) = A(i,k)/A(k,k);
        b(i) = b(i) - A(i,k)*b(k);
        for j = k+1:n
            A(i,j) = A(i,j) - A(i,k)*A(k,j);
        end
    end
end
        
% Extrahieren der Matrizen L und U
U = zeros(n,n);
for k=1:n
    for i=1:k
        U(i,k) = A(i,k); % oberer Dreiecksanteil von A
    end
end
L = eye(n);
for k=1:n
    for i=k+1:n
        L(i,k) = A(i,k); % (strikter) unterer Dreiecksanteil von A
    end
end

y = b;

x = zeros(n,1);
for i = n:-1:1 % Rueckwaertssubst.
    sum = 0;
    for j = i+1:n
        sum = sum+U(i,j)*x(j);
    end
    x(i) = (y(i) - sum)/U(i,i);
end

fprintf('Residuum fuer Matrixfakt. = %g\n',norm(L*U - A_copy,'fro'));

