% script file: gFHSVDDemo.m

% Demonstrates general Hankel Lanczos bidiagonalization and SVD.

%

% Dependency

%    ./gFHLanMPO.m   general fast hankel Lanczos bidiagonlization with MPO 

%    ./gFHLanMPOR.m  general fast hankel Lanczos bidiagonlization with MPO 

%                    and reset

%    ./gLanPO.m      general Lanczos with PO

%    ./gLanMPO.m     general Lanczos bidiagonalization with MPO

%    ./gLanMPOR.m    general Lanczos bidiagonalization with MPO and reset

%    ./appsvd.m      symmetric SVD of a complex-symmetric and bidiagonal matrix

%

% K. Browne       McMaster Univ.  March 2008



% input the row and column dimensions

m = input('Enter the row dimension: m = ');

n = input('Enter the column dimension: n (n<=m) = ');

while (n > m)

    disp('Error: column dim is greater than row dim.')

    n = input('Enter another column dimension: ');

end





% generate the first col and last row of Hankel

col = (2*rand(m,1) - ones(m,1)) + j*(2*rand(m,1) - ones(m,1));

row = (2*rand(n,1) - ones(n,1)) + j*(2*rand(n,1) - ones(n,1));

row(1) = col(m);	% last element of col wins over

			% first element of row

A = hankel(col,row);





% general Lanczos bidiagonalization with partial orthogonalization.

[a,b,V1,U1] = gLanPO(A,ones(n,1));



% general Lanczos bidiagonalization with modified partial orthogonalization

% [a,b,V1,U1] = gLanMPO(A,ones(n,1));



% general Lanczos bidiagonalization with modified partial orthogonalization 

% and restart.

% [a,b,V1,U1] = gLanMPOR(A,ones(n,1));



% general fast Lanczos bidiagonalization of Hankel using modified

% partial orthogonalization 

% [a,b,V1,U1] = gFHLanMPO(col,row,ones(n,1));



% general fast Lanczos bidiagonalization of Hankel using modified

% partial orthogonalization with restart

% [a,b,V1,U1] = gFHLanMPOR(col,row,ones(n,1));



% report errors

B = (diag(a) + diag(b,1));

tmp = norm(eye(n) - U1'*U1);

fprintf('\nError in orthogonality of U in bidiagonalization: %E', tmp);

tmp = norm(eye(n) - V1'*V1);

fprintf('\nError in orthogonality of V in bidiagonalization: %E', tmp);

tmp = norm(A - U1*B*V1');

fprintf('\nError in bidiagonalization: %E', tmp);





% perform bidiagonal SVD

fnrm = hankelfnrm(col,row);   % F-norm of Hankel

tolsd = sqrt(eps)*fnrm/(m*n);

                  % tolerance for small superdiagonal elements

[s,U2,V2,smlflag] = appsvd(a,b,tolsd);





% report errors comparing with builtin svd

sv = svd(A);

tmp = norm(sv - s);

fprintf('\nError in singular values: %E', tmp);



% if we have U and V, compute all the errors

if (~smlflag)

    tmp = norm(B - (U2 * diag(s) * V2'));

    fprintf('\nError in bidiagonal SVD: %E', tmp);

    V = V1 * V2;

    U = U1 * U2;

    tmp = norm(A - (U * diag(s) * V'));

    fprintf('\nError in combined bidiagonlization and SVD: %E\n', tmp);

end