% --- Icke-nolljägare (med illustrativa möjligheter) ---
% Chalmers Tekniska Högskola, 2002-10-09
% Numerisk Linär Algebra
%
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
    % nollar xj i 'i' kolonnen
    xi = A(i+1,i);
    xj = A(i+2,i);

    C = xi/sqrt(xi^2+xj^2);                 % cosinustermen
    S = -xj/sqrt(xi^2+xj^2);                % sinustermen
    R = eye(m);
    R(i+1:i+2,i+1:i+2) = [C -S; S C];       % Givens-rotationsmatrisen

    spy(A)                                  % Visar hur R-matrisen ser
    hold on                                 % ut i jämförelse med
    spy(R,'r*')                             % A-matrisen. Bör endast
    hold off                                % användas då
    drawnow                                 % A = round(R*A*R')
    
    A = round(R*A*R');                      % för en exakt lösning, ändra till A = R*A'R
    j = i;
    
    % jagar bort skräptermen
    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)                              % Visar hur R-matrisen ser
        hold on                             % ut i jämförelse med
        spy(R,'r*')                         % A-matrisen. Bör endast
        hold off                            % användas då
        drawnow                             % A = round(R*A*R')
    
        A = round(R*A*R');                  % för en exakt lösning, ändra till A = R*A'R


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

sort(eig(full(A)))-sort(eig(full(Aorg)))    % validerar egenvärdens bevarande

% av Markus Burak & Tobias Tivert