function sigit = mysvdplot ( nit, A )
% mySVDPlot
%
% Computation of Singular Values, Plot version
%
% Parameters
%   nit    (input)  : number of iteration steps in QR algorithm
%   A      (input)  : matrix A with m rows, n columns, m>=n
%   sigit  (input)  : diagonal entries of the bidiagonal iterates in SVD
%
% Example
%   nit = 60;  A = rand ( 36, 24 );
%   sigit = mysvdplot ( nit, A );
%
% Reference
%   P. Deuflhard, A. Hohmann: Numerische Mathematik 1, 3. Auflage.
%      Walter de Gruyter Berlin New York, 2002, Abschnitt 5.4.
%
% Author :      Martin Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Nov 24, 2009

% -> parameters of graphical output
CurFontSize  = 14;
CurLineWidth =  2;

cmin = -3;
cmax = log10 ( norm ( A ) );

clf

% -> initialization of output data
[ m, n ] = size ( A );
if m<n
    error ( 'Matrix A (input): Number n of columns must not exceed number m of rows' )
end

p = min ( m, n );
sigit = zeros ( nit+1, p );

% -> initialization of extended spy plot
root = 0;
screensize = get ( root, 'ScreenSize' );
left   = screensize(1);
bottom = screensize(2);
width  = screensize(3);
height = screensize(4);

if isempty (get(root,'Children')),
  h = figure;
else
  h = gcf;
  clf;
end;
set ( h, 'Position', [ left+0.15*width bottom+0.3*height 0.8*width 0.6*height ] );

[ UU, SS, VV ] = svd ( A );

X = zeros ( 4, m*n );
Y = zeros ( 4, m*n );
C = zeros ( 1, m*n );

for ic=1:m
    for jc=1:n
        icol = n*(ic-1)+jc;
        X(:,icol) = jc + 0.5*[ -1 -1 1  1 ]';
        Y(:,icol) = ic + 0.5*[ -1  1 1 -1 ]';
        C(icol)   = log10 ( eps + abs ( A(ic,jc) ) );
    end
end

subplot ( 1, 2, 1 )
matpic = fill ( X, Y, C );
axis ( [ 0 n+1 0 m+1 ] )
set ( gca, 'FontSize', CurFontSize );
set ( gca, 'CLim', [ cmin cmax ] );
set ( gca, 'YDir'              , 'reverse' )
set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] )
set ( gca, 'DataAspectRatio'   , [ 1 1 1 ] )
climv = get ( gca, 'CLim' );
cb = colorbar;
set ( cb, 'FontSize', CurFontSize )
set ( cb, 'YTick', climv )
set ( cb, 'YTickLabel', { ' 0', ' max' } )
cm = colormap;
drawnow

% -> bi-diagonalization
for k=1:p
    
    % -> Householder reflection from the left
    if A(k,k)>=0
        sa = - norm ( A(k:end,k) );
    else
        sa = norm ( A(k:end,k) );
    end
    fac = sa * ( sa - A(k,k) );
    v = [ zeros(k-1,1); A(k,k)-sa; A((k+1):end,k) ];
    % A(k,k) = sa;
    for j=k:n
        A(k:end,j) = A(k:end,j) - v(k:end)'*A(k:end,j) * v(k:end) / fac;
    end

    for ic=1:m
        for jc=1:n
            set ( matpic(n*(ic-1)+jc), 'CData', log10(eps+abs(A(ic,jc))) );
        end
    end
    drawnow
    
    % -> Householder reflection from the right
    if k<p
        if A(k,k+1)>=0
            sa = - norm ( A(k,(k+1):end) );
        else
            sa = norm ( A(k,(k+1):end) );
        end
        fac = sa * ( sa - A(k,k+1) );
        v = [ zeros(1,k), A(k,k+1)-sa, A(k,(k+2):end) ]';
        % A(k,k+1) = sa;
        for i=k:m
            A(i,(k+1):end) = A(i,(k+1):end) - A(i,(k+1):end)*v((k+1):end) * v((k+1):end)' / fac;
        end
    end

    % -> update extended spy plot
    for ic=1:m
        for jc=1:n
            set ( matpic(n*(ic-1)+jc), 'CData', log10(eps+abs(A(ic,jc))) );
        end
    end
    drawnow

end

it = 0;
sigit(it+1,:) = abs ( diag ( A ) );

% -> initialize error plot 
subplot ( 1, 2, 2 )
svpic = semilogy ( 0:nit, NaN*ones(size(sigit)), 'LineWidth', CurLineWidth );
set ( gca, 'FontSize', CurFontSize );
for ic=1:p
    ydat = get ( svpic(ic), 'YData' );
    ydat(1:(it+1)) = abs(sigit(1:(it+1),ic)-SS(ic,ic));
    set ( svpic(ic), 'YData', ydat )
    colic = ( log10(eps+abs(A(ic,ic))) - cmin ) / ( cmax - cmin );
    icm = max ( 1, min ( length(cm(:,1)), 1+floor(colic*(length(cm(:,1))-1)) ) );
    set ( svpic(ic), 'Color', cm(icm,:) )
end  
axis ( [ 0 nit 1.0e-15 1.0 ] )
xlabel ( 'k' )
ylabel ( '| \sigma_i - B_{ii}^{(k)} |' )
title ( 'SVD by QR-algorithm' )
drawnow

% -> QR algorithm
for it=1:nit
    
    if abs(A(1,1))>abs(A(1,2))
        tau = A(1,2) / A(1,1);
        c   = 1 / sqrt ( 1 + tau^2 );
        s   = -c * tau;
    else
        tau = A(1,1) / A(1,2);
        s   = 1 / sqrt ( 1 + tau^2 );
        c   = -s * tau;
    end
    
    A(1:2,1:2) = [ A(1,1:2); 0 A(2,2) ] * [ c s ; -s c ];
    for ic=1:2
        for jc=1:2
            set ( matpic(n*(ic-1)+jc), 'CData', log10(eps+abs(A(ic,jc))) );
        end
    end
    drawnow
    
    for i=1:(p-1)

        if abs(A(i,i))>abs(A(i+1,i))
            tau = A(i+1,i) / A(i,i);
            c   = 1 / sqrt ( 1 + tau^2 );
            s   = -c * tau;
        else
            tau = A(i,i) / A(i+1,i);
            s   = 1 / sqrt ( 1 + tau^2 );
            c   = -s * tau;
        end
    
        if i<(p-1)
            A(i+(0:1),i+(0:2)) = [ c -s ; s c ] * [ A(i,i+(0:1)) 0 ; A(i+1,i+(0:2)) ];
            for ic=i+(0:1)
                for jc=i+(0:2)
                    set ( matpic(n*(ic-1)+jc), 'CData', log10(eps+abs(A(ic,jc))) );
                end
            end
            drawnow
            if abs(A(i,i+1))>abs(A(i,i+2))
                tau = A(i,i+2) / A(i,i+1);
                c   = 1 / sqrt ( 1 + tau^2 );
                s   = -c * tau;
            else
                tau = A(i,i+1) / A(i,i+2);
                s   = 1 / sqrt ( 1 + tau^2 );
                c   = -s * tau;
            end
            A(i+(0:2),i+(1:2)) = [ A(i+(0:1),i+(1:2)); 0 A(i+2,i+2) ] * [ c s ; -s c ];
            for ic=i+(0:2)
                for jc=i+(1:2)
                    set ( matpic(n*(ic-1)+jc), 'CData', log10(eps+abs(A(ic,jc))) );
                end
            end
            drawnow
        else
            A(i+(0:1),i+(0:1)) = [ c -s ; s c ] * [ A(i,i+(0:1)); A(i+1,i+(0:1)) ];
            for ic=i+(0:1)
                for jc=i+(0:1)
                    set ( matpic(n*(ic-1)+jc), 'CData', log10(eps+abs(A(ic,jc))) );
                end
            end
            drawnow
        end
    end
    
    sigit(it+1,:) = abs ( diag ( A ) );

    % -> update error plot
    for ic=1:p
        ydat = get ( svpic(ic), 'YData' );
        ydat(1:(it+1)) = abs(sigit(1:(it+1),ic)-SS(ic,ic));
        set ( svpic(ic), 'YData', ydat )
        colic = ( log10(eps+abs(A(ic,ic))) - cmin ) / ( cmax - cmin );
        icm = max ( 1, min ( length(cm(:,1)), 1+floor(colic*(length(cm(:,1))-1)) ) );
        set ( svpic(ic), 'Color', cm(icm,:) )
    end
    drawnow

end



