%% Komplexní vlastní čísla
%% QR algoritmus bez posunutí

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

% Další 2 iterace

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

%
%% Rayleighovo posunutí

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

% Další 4 iterace

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

%
% Submatice:

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

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

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

% Další 2 iterace

A2=QRAW(A2,1)
A2=QRAW(A2,1)

%

A3=A2(1:2,1:2)
A4=QRAW(A3,1)
A4=QRAW(A3,1)
eig(A4)

%

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

