
%    This script tests the randomized low rank SVD code randSVD
%
     m = 500; n = 50;

%
%    Generate a random orthogonal matrix Uo of order m by n , m >> n
%
     [Uo,R] = qr(randn(m,n));
%
%    Do full SVD to test subspace convergence
%
     [V,R] = qr(randn(n)); % generate random orthogonal n by n  V
%
%    Set left singular vectors U to leading n cols. of Uo
%
     U = Uo(:,1:n);

     s = ones(n,1);
     t = 2;
     alpha = -1;
     beta = n;
%
%    Generate a vector s of singularvalues
%
%        case test == 'lin'  sets s(j) = alpha*j + beta
%
%        case test == 'exp'  sets s(j+1) = 1/(2^j) j = 0:1:n-1 

     test = input('  lin or exp ')
     for j = 2:n, s(j) = 1/t; t = t*2; end

     if (test == 'lin'), 

         for j = n:-1:1, s(j) = alpha*j + beta; end; 

         s = s/s(1); % scale singular values to be in [0 , 1]
     end

     A = U*diag(s)*V';

     k = input(' value of k = ')
     q = input(' value of q = ')

     [U,S,V] = randSVD(A,k,q,true);  % Iterates with AA'

%    [U,S,V] = randSVDApA(A,k,q,true);    % Iterates with A'A

%
     disp(' SVD approx error vs s(k+1) ')

     Dnorm_Sigma = [norm(A - U*S*V') s(k+1)]
     ss = diag(S);

%
%    compute distance between approximate dominant subspace Range(U(:,1:k))
%    and the true subspace  Range(Uo(:,1:k))  (see G and V pp76-77 )
%

     disp(' Distance between Range(U(:,1:k)) and Range (Uo(:,1:k))')

     SSdistance = norm(U(:,1:k)'*Uo(:,k+1:m))

     semilogy(s,'ro', 'linewidth',1.5) 
     hold on
     semilogy(ss,'.', 'linewidth',1.5) 
     hold off

     title(' Compare Approx vs True Singular Values', 'fontsize',24)
     xlabel(' Index ','fontsize', 18)
     ylabel(' Singular Values ','fontsize', 18)
     legend('true svals', 'approx svals')



     disp('Approx vs true singular values')

     disp ( '    approx        true         diff ')
%    Sdiff = [ss s(1:length(ss)) abs(ss - s(1:length(ss)))]
     disp([ss s(1:length(ss)) abs(ss - s(1:length(ss)))])
