%% QR algoritmus bez posunutí

format short
A=[11 -6 4 -2; 4 1 0 0; -9 9 -6 5; -6 6 -6 7]
eig(A)
AA=hess(A);
for k=1:10
    [Q,R]=qr(AA);
    AA=R*Q;
end
AA
AA(4,3)

%
%% Rayleighovo posunutí

format short
A=[11 -6 4 -2; 4 1 0 0; -9 9 -6 5; -6 6 -6 7];
AA=hess(A);
I=eye(size(AA));
for k=1:10
    mi=AA(end,end);
    MI=mi*I;
    [Q,R]=qr(AA-MI);
    AA=R*Q+MI;
end
AA
eig(AA(2:3,2:3))

%
% Další krok:

mi=AA(end,end);
MI=mi*I;
[Q,R]=qr(AA-MI);
AA=R*Q+MI
eig(AA(2:3,2:3))

%
% Submatice:

AA=AA(1:3,1:3)
I=eye(size(AA));
for k=1:10
    mi=AA(end,end);
    MI=mi*I;
    [Q,R]=qr(AA-MI);
    AA=R*Q+MI;
end
AA
AA(2,1)

%
%% Wilkinsonovo posunutí
% pro předchozí matici dává obdobné výsledky jako Rayleghovo posunutí

A1=hess(A);
A2=QRAW(A1,10)

%

A3=A2(1:3,1:3);
A4=QRAW(A3,10)

%
% Příklad, kdy Rayleighovo posunutí nefunguje

A=[0 1;1 0];
eig(A)
AA=A;

% 1 krok QR algoritmu - matice zůstává beze změny

[Q,R]=qr(AA);
AA=R*Q

%
% Rayleighovo posunutí nelze použít, zkusíme jinou hodnotu
% 1 iterace

AA=A;
mi=0.1;
I=eye(size(AA));
MI=mi*I;
[Q,R]=qr(AA-MI);
AA=R*Q+MI

%

for k=1:5
    mi=AA(end,end);
    MI=mi*I;
    [Q,R]=qr(AA-MI);
    AA=R*Q+MI
end
%AA

%
% Wilkinsonovo posunutí

A1=A;
A2=QRAW(A1,1)
A2(2,1)

%

function [AA] = QRAW(A,k)
% function [AA] = QRAW(A,k)
% k kroku QR algoritmu s Wilkinsonovym posunutim   
%
n=length(A);
I=eye(n);
for ii=1:k
    AA=A(n-1:n,n-1:n);
    d=det(AA);
    t=trace(AA);
    if t^2<4*d      % komplexni vlastni cisla
        M=A*A-t*A+d*I;
        [Q,R]=qr(M);
        A=Q'*A*Q;
    else
        a11=AA(1,1);
        a12=AA(1,2);
        a21=AA(2,1);
        a22=AA(2,2);
        delta=(a11-a22)/2;
        if delta==0 
            sd=1;
        else
            sd=sign(delta);
        end
        mi=a22-sd*a12*a21/(abs(delta)+sqrt(delta^2+a12*a21));
        [Q,R]=qr(A-mi*I);
        A=R*Q+mi*I;
    end
end
AA=A;
end

