% --- Bulge chasing (with possible illustrations) ---
% Chalmers University of Technology, 2002-10-09
% Numerical Linear Algebra
% By Markus Burak and Tobias Tivert
clear all
clc

m=40;
e = ones(m,1);
A = 1000*full(spdiags([1*e 2*e 4*e 2*e 1*e], [-2:2], m, m));
n = length(A);
Aorg = A;

for i=1:n-2
    % zero-out xj in i'th column
    xi = A(i+1,i);
    xj = A(i+2,i);

    C = xi/sqrt(xi^2+xj^2);                 % cos-term
    S = -xj/sqrt(xi^2+xj^2);                % sin-term
    R = eye(m);
    R(i+1:i+2,i+1:i+2) = [C -S; S C];       % Givens' rotation matrix

    spy(A)                                  % Shows how the R-matrix looks
    hold on                                 % compared with the
    spy(R,'r*')                             % A-matrix.
    hold off                                % To be used only if
    drawnow                                 % A = round(R*A*R')
    
    A = round(R*A*R');                      % for exact solution - change to A = R*A'R
    j = i;
    
    % bulge chasing
    while j < n-3,

        xi = A(j+3,j+1);
        xj = A(j+4,j+1);
        C = xi/sqrt(xi^2+xj^2);
        S = -xj/sqrt(xi^2+xj^2);
        R = eye(m);
        R(j+3:j+4,j+3:j+4) = [C -S; S C];

        spy(A)                              % shows how the R-matriix looks
        hold on                             % compared to the
        spy(R,'r*')                         % A-matrix.
        hold off                            % iTo be used only if
        drawnow                             % A = round(R*A*R')
    
        A = round(R*A*R');                  % for exact solution - change to A = R*A'R


        j = j + 2;
    end
end
spy(A)                                      
drawnow                                     

sort(eig(full(A)))-sort(eig(full(Aorg)))    % check the preservation of eigenvalues
