%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Aufgabe 3:
%   Schreiben Sie eine Funktion solv, der man eine Matrix A und einen 
%   Vektor b übergibt. 
%    (1) Bestimmen Sie eine LR-Zerlegung PA=LR und lösen Sie die 
%        Gleichungssysteme Lu=Pb und Ry=u ohne die Matrizen L und R zu 
%        invertieren.
%    (2) Vergleichen Sie den Vektor y mit der Lösung x von Ax=b, die sie 
%        über den Matlab \ Operator berechnen.
%    (3) Geben Sie die Lösung x und den Fehlervektor |y-x|, sowie den 
%        durchschnittlichen Fehler mean(|y-x|) und die 
%        1-norm-Konditionszahl von A in dieser Reihenfolge zurück.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,fv,df,k]=BSP_3(A,b)
% teste die Funktion mit: A=rand(3,3);b=rand(3,1);[x,fv,df,k]=BSP_3(A,b);

m=length(b);

%(1.1)
[L,R,P]=lu(A); % LR-Zerlegung PA=LR => LRy=Pb
%(1.2) berechene u=L\(P*b) mittels Vorwärtssubstitution
Pb=P*b;
u=zeros(m,1);
for i=1:m
    %u(i):=(Pb(i)-L(i,J)*u(J))/L(i,i), J={1,...,i-1}
    J=1:i-1;
    u(i)=(Pb(i)-L(i,J)*u(J))/L(i,i); 
end
%TEST: L*u-Pb sollte 0 sein 
%(1.2) berechene y=R\u; mittels Rückwärtssubstitution
y=zeros(m,1);
for i=m:-1:1
    %y(i):=(u(i)-R(i,J)*y(J))/R(i,i), J={m,...,i+1}
    J=m:-1:i+1;
    y(i)=(u(i)-R(i,J)*y(J))/R(i,i);
end
%TEST: R*y-u und A*y-b sollte 0 sein 
disp(['(1) Lösung y=' mat2str(y)]);

%(2)
x=A\b;  % lösung x von Ax=b, mittels den Matlab \ Operator
disp(['(2) Lösung x=' mat2str(x)]);

%(3)
fv=abs(y-x); % Fehlervektor
df=mean(fv); % durchschnittlicher Fehler
k=cond(A,1); % 1-norm-Konditionszahl von A

disp('(2) Lösung und Rückgabe');
disp(['    fv=' mat2str(fv)]);
disp(['    df=' mat2str(df)]);
disp(['    k='  mat2str(k)]);

