%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%                      EXPLICITNI QR-ALGORITMUS: 
%               konvergence, preprocesing, Rayleighuv shift
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

A=gallery(3);                        % testovaci matice s vl. c. 1, 2, 3
eig(A)
 
%% naprogramujte zakladni explicitni QR-algoritmus
 
[R,err,lambda_n]=explQR(A,50);       % pozorujte matici A v prubehu iteraci

subplot(2,1,1)
semilogy(err)                        % pozorujte velikost prvku A(n,n-1)
                                     % -> konverguje, diverguje, konverguje
subplot(2,1,2)
plot(lambda_n,'*')                   % aproximace nejmensiho vl.c.                                   
                                     % -> vykresele bez prvnich 5-ti iteraci
                                     % -> nejprve konv. k jinemu vl. cislu!
                                     
%% pridejte preprocesing 

%help hess                           % podobnostni transformace na horni
[P,H] = hess(A);                     % Hessenberguv tvar
%figure(2), spy(H)

[R,err,lambda_n]=explQR(H,50);       % QR alg. na horni Hes. matici             

figure(3)
subplot(2,1,1)
semilogy(err)  
subplot(2,1,2)
plot(lambda_n,'*')                   % konverguje k vl. cislu 1    

% sledujte konvergenci ostatnich poddiagonalnich prvku k 0 a tim padem 
% diagonalnich prvku k vl c. A


%% pridejte do explicitniho QR-algoritmu Rayleighuv shift

[R1,err1,lambda1]=explQR_shift(H,50);     

subplot(2,1,1)
hold on
semilogy(err1,'m')  
subplot(2,1,2)
hold on
plot(lambda1,'mo')            % konverguje mnohem rychleji k nejm. vl.c.

% sledujte konvergenci ostatnich poddiagonalnich prvku k 0 a tim padem 
% diagonalnich prvku k vl c. A

%% pridejte zastavovaci kriterium, redukci rozmeru a rekurzi

%% testujte dalsi matice:

% Priklad 2: D = diag([4 3 2 1]); S=rand(4); A = S*D/S;
% -> konverguje vse, se shiftem nejrychleji
% -> shifted QR bez inicializacnich iteraci konverguje k jinemu poradi vl.c.

% Priklad 3: D = diag([4 3 -1 1]); S=rand(4); A = S*D/S; 
% -> zakladni alg. diverguje, alg. se shiftem konverguje  

% Priklad 4: A=randn(50); A = A+A';
% -> symetricka matice
% -> H je tridiagonalni, R konverguje k diagonalni ... surf(R)