function Keeling

% make a matrix A
    m = 10;
    n = 3;
    A = rand(m,n);
% transform A to bidiagonal form B so that P*A*Q = (B,0)^T
    [P,B,Q] = ToBiDiag(A);
% check the result
    PAQ = P*A*Q;
    PAQ = PAQ(1:n,1:n);
    disp(['residual in bidiagonal transformation = ', ...
        num2str(norm(PAQ-B))]);
% compute singular values of B
    S = BiDiagQR(B);
% check the result
    disp(['residual in singular value determination = ', ...
        num2str(norm(sort(S)-sort(svd(A))))]);

end

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

function [P,B,Q] = ToBiDiag(A);

    [m,n] = size(A);
    if (m<n)
        error('Repeat with A'' so that m>n.');
    end

    d = A(:,1);
    C = A(:,2:n);
    B = zeros(n,n);

    P = eye(m);
    Q = eye(n);
    B = zeros(n,n);

    for k=2:(n+1)
        [u,theta,sigma] = Householder(d);
        B(k-1,k-1) = -sigma;
        if (k>2)
            B(k-2,k-1) = a;
        end
        for l=1:(n-k+1)
            sm = dot(C(:,l),u);
            C(:,l)  = C(:,l) - sm*u/theta;
        end
        for l=1:m
            sm = dot(P((k-1):m,l),u);
            P((k-1):m,l) = P((k-1):m,l) - sm*u/theta;
        end
        if (k < n)
            c = C(1,:);
            [u,theta,tau] = Householder(c);
            a = -tau;
            for l=2:(m-k+2)
                sm = dot(C(l,:),u);
                C(l,:)  = C(l,:) - sm*u/theta;
            end
            for l=1:n
                sm = dot(Q(l,k:n),u);
                Q(l,k:n) = Q(l,k:n) - sm*u/theta;
            end
        end
        if (k == n)
            a = C(1,1);
        end
        if (k < (n+1))
            d = C(2:(m-k+2),1);
            C = C(2:(m-k+2),2:(n-k+1));
        end
    end
end

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

function S = BiDiagQR(B);

    n = size(B,1);
    kmax = 1000;
    tol  = 1.0e-10;

    x = zeros(n,1); g=x; y = zeros(n-1,1); c=y; s=y; d=y; e=zeros(n-2,1);

    for k=1:kmax

        BB = B'*B;
        a  = diag(BB);
        b  = diag(BB,-1);

        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);
        R = diag(g) + diag(d,1) + diag(e,2);
        Q = eye(n);
        for j=1:(n-1)
            Q1 =  Q(:,j)*c(j) + Q(:,j+1)*s(j);
            Q2 = -Q(:,j)*s(j) + Q(:,j+1)*c(j);
            Q(:,j  ) = Q1;
            Q(:,j+1) = Q2;
        end

        BQ = B*Q;
        a = diag(BQ);
        b = diag(BQ,-1);
        x(1) = a(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) = s(j)*a(j+1);
            x(j+1) = c(j)*a(j+1);
        end
        g(n) = x(n);
        B = diag(g) + diag(d,1);

        if (norm(d) < tol*norm(g))
            break;
        end
    end
    S = abs(g);
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

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