function Keeling

% make a matrix A
    n = 10;
    A = rand(n,n);
    A = A'*A + eye(n);
% transform A to tridiagonal form H so that A = Q*H*Q'
    [H,Q] = ToTriDiag(A);
% check the result
    disp(['residual in tridiagonal transformation = ', ...
        num2str(norm(A-Q*H*Q'))]);
% compute eigenvalues of H
    E = TriDiagQR(H);
% check the result
    disp(['residual in eigenvalue determination = ', ...
        num2str(norm(sort(E)-sort(eig(A))))]);

end

%%--------------------------------------------------------

function [H,Q] = ToTriDiag(A)
    n = size(A,1);
    H = A;
    Q = eye(n);
    for k=1:(n-2)
        [u,theta,sigma] = Householder(H((k+1):n,k));
        for l=k:n
            sm = dot(H((k+1):n,l),u);
            H((k+1):n,l) = H((k+1):n,l) - sm*u/theta;
        end
        for l=1:n
            sm = dot(Q(l,(k+1):n),u);
            Q(l,(k+1):n) = Q(l,(k+1):n) - sm*u'/theta;
        end
        for l=k:n
            sm = dot(H(l,(k+1):n),u);
            H(l,(k+1):n) = H(l,(k+1):n) - sm*u'/theta;
        end
    end
end

%%--------------------------------------------------------

function a = TriDiagQR(H);
    n = size(H,1);
    a = diag(H);
    b = diag(H,1);
    tol  = 1.0e-10;
    kmax = 1000;
    for k=1:kmax
        x(1) = a(1);
        y(1) = b(1);
        for j=1:(n-1)
            g(j) = sqrt(x(j)^2 + b(j)^2);
            c(j) = x(j)/g(j);
            s(j) = b(j)/g(j);
            d(j) = c(j)*y(j) + s(j)*a(j+1);
            x(j+1) = -s(j)*y(j) + c(j)*a(j+1);
            if (j ~= (n-1))
                e(j) = s(j)*b(j+1);
                y(j+1) = c(j)*b(j+1);
            end
        end
        g(n) = x(n);
        a(1) = s(1)*d(1) + c(1)*g(1);
        b(1) = s(1)*g(2);
        for j=2:(n-1)
            a(j) = s(j)*d(j) + c(j)*c(j-1)*g(j);
            b(j) = s(j)*g(j+1);
        end
        a(n) = c(n-1)*g(n);
        if (norm(b) < tol)
            break;
        end
    end
end

%%--------------------------------------------------------

function [u,theta,sigma] = Householder(x)
    s = @(t) double(t>=0) - double(t<0);
    xinf  = norm(x,inf);
    v     = x/xinf;
    rho   = s(v(1))*norm(v,2);
    u     = v;
    u(1)  = u(1) + rho;
    theta = rho*u(1);
    sigma = rho*xinf;
end

%%--------------------------------------------------------
