function Blatt6
format long
% Aufgabe 6.4 - 1. Beispiel
A1 = [2,0,1;6,-1,6;-4,-5,15];
b1 = [5;10;15];
% LR Zerlegung
[L1,R1] = LR_decomp(A1);
% Ergebnis mit Vorwärts/Rückwärtsauflösen
x1 = Vorwaertsaufloesen(L1,b1);
Erg1 = Rueckwaertsaufloesen(R1,x1);
disp(Erg1)
% gegentesten des Vorwärts/Rückwärtsauflösens
x1test = L1\b1;
Erg1_test = R1\x1test;
disp(Erg1_test)
% gegetesten des gesamten Ergebnisses
Matlab_Erg1 = A1\b1;
disp(Matlab_Erg1)

% 2. Beispiel
A2 = [2,3,0;2,3,1;2,1,1];
b2 = [1;1;1];

[L2,R2] = LR_decomp(A2);
% man sieht, dass für j = 2 das Pivotelement = 0 ist. Daher ist es wichtig die LR-Zerlegung mit Pivotisierung zu implementieren!
end

function [L,R] = LR_decomp(A)
  % LR-Zerlegung ohne Pivotsuche nach Algorithmus 1 der Vorlesung
  n = size(A,1); % A ist quadratisch  
  % Rest entspricht dem Pseudocode der VL
  for j = 1:n-1
    if abs(A(j,j)) < 1e-12 % Falls das Diagonalelement aufgrund von Rechenungenauigkeit numerisch nicht exakt 0 ist, führt das zu Problemen. So fängt man dieses Problem ab
      error('Pivotelemnt gleich Null!')
    else
      for k = j+1:n
        A(k,j)= A(k,j)/A(j,j); % l_kj
        for i = j+1:n
          A(k,i) = A(k,i) - A(k,j)*A(j,i); %A_ki^j+1
        end
%        A(k,j+1:end) = A(k,j+1:end) - A(k,j)*A(j,j+1:end); % so könnte man sich die 3. for-Schleife noch sparen...
      end
    end
  end
  L = tril(A,-1)+eye(n,n);
  R = triu(A);
end

function x1 = Vorwaertsaufloesen(L,b)
  % Vorwärtsauflösen
  n = size(L,1);
  x1 = zeros(n,1);
  for j = 1:n
    Abzug = sum(L(j,1:j-1)*x1(1:j-1)); % die Elemente in L links von der Diagonalen müssen gewichtet mit den bisherigen x abgezogen werden
    x1(j) = b(j) - Abzug; % keine Skalierung der b(j), da Einsen auf der Diagonale sind
  end
end

function [x] = Rueckwaertsaufloesen(R,x1)
  % Rückwärtsauflösen - initialisiere auch x, da eben Rückwärtsaufgelöst wird
  n = size(R,1);
  x = zeros(n,1);
  for j = 1:n
    Abzug = sum(R(end-j+1,end-j+2:end)*x(end-j+2:end)); % Die Elemente von R rechts der Diagonalen müssen abgezogen werden
    x(end-j+1) = (x1(end-j+1)-Abzug)/R(end-j+1,end-j+1); % skalierung mit dem Diagonalelement
  end
end