% Skript zur Musterlösung
%
% Numerische Mathematik, SS 2013
%
% Fehlerschranken mit dem Hager-Algorithmus
%
% Ingolf Busch, Juni 2013

% Initialisierung
n=5:100;
err_pp=zeros(size(n));
err_pp_1=zeros(size(n));
err_pp_2=zeros(size(n));

err_cp=zeros(size(n));
err_cp_1=zeros(size(n));
err_cp_2=zeros(size(n));

% Funktionen für Matrix-Vektor-Produkt

mv  = @(A,x) my_solve_lr(A,x,false);
mvt = @(A,x) my_solve_lr(A,x,true);

mvg  = @(A,x,g) my_solve_lr(A,g.*x,false);
mvgt = @(A,x,g) g.*my_solve_lr(A,x,true);


for i = 1:length(n)
    
    %Matrix generieren
    c=10^(14/(n(i)-1));
    
    D=diag(c.^(0:(n(i)-1)));
    B = eye(n(i)) + randn(n(i))*1e-7;
    
    A = D*B;
    
    % rechte Seite erzeugen
    x = randn(size(A,1),1);
    b = A*x;
    
    % Spaltenpivotsuche
      
    S=my_lr_pp(A);
    
    x_pp = my_solve_lr(S,b,false);
    err_pp(i)=norm(x-x_pp,Inf)/norm(x_pp,Inf);
    r_pp = abs(b-A*x_pp);
    
    err_pp_1(i) = my_hager(S,size(A,1),mvt,mv,100)*norm(r_pp,Inf)/norm(x_pp,Inf);
    err_pp_2(i) = my_hager(S,size(A,1),@(A,x) mvgt(A,x,r_pp),@(A,x) mvg(A,x,r_pp),100)/norm(x_pp,Inf);
    
    % vollständige Pivotsuche
    
    S=my_lr_cp(A);
    x_cp = my_solve_lr(S,b,false);
    err_cp(i)=norm(x-x_cp,Inf)/norm(x,Inf);
    r_cp=abs(b-A*x_cp);
    
    err_cp_1(i) = my_hager(S,size(A,1),mvt,mv,100)*norm(r_cp,Inf)/norm(x_cp,Inf);
    err_cp_2(i) = my_hager(S,size(A,1),@(A,x) mvgt(A,x,r_cp),@(A,x) mvg(A,x,r_cp),100)/norm(x_cp,Inf);    
end

figure;
loglog(10.^[-18 -2],10.^[-18 -2],'b-',10.^[-18 -2],(10.^(1:15))'*(10.^[-18 -2]),'b--');
hold on;
h=loglog(err_pp,err_pp_1,'go',err_cp,err_cp_1,'r+');
title('Fehlerschätzer 1');
xlabel('wahrer Fehler');
ylabel('geschätzter Fehler');
legend(h,'partial pivot','complete pivot','Location','SouthEast');
axis(10.^[-18 -2 -18 -2]);

figure;
loglog(10.^[-18 -2],10.^[-18 -2],'b-',10.^[-18 -2],(10.^(1:15))'*(10.^[-18 -2]),'b--');
hold on;
h=loglog(err_pp,err_pp_2,'go',err_cp,err_cp_2,'r+');
title('Fehlerschätzer 2');
xlabel('wahrer Fehler');
ylabel('geschätzter Fehler');
legend(h,'partial pivot','complete pivot','Location','SouthEast');
axis(10.^[-18 -2 -18 -2]);