


function [U,S,V] = randSVD(A,k,q,orth);
%
%   This function computes an approximate truncated
%   SVD  for  A ~ Uk Sk Vk' of rank 2k.  The accuracy
%   should be 
%
%             || A - Uk Sk Vk' || = O(sigma_k+1) 
%   
%             with sigma_k+1 the first neglected singular value
%
%   Uses Y = (A A')^q A W   to get approx basis Y for Range(A)
%
%
%   Input:  
%          A      an m by n matrix, m .ge. n 
%
%          k      a positive integer (k << n)
%
%          q      the number of power iterations with AA'
%
%          orth   a logical variable.  Value .true.  indicates
%                 re-orthogonalization after each power step,
%                 while .false. indicates no re-orthogonalization
%
%   Output:
%          U  an m by 2k  orthogonal matrix
%          
%          S  a  2k by 2k diagonal matrix of 
%             approximate singular values
%
%          V  a 2k by n  orthogonal matrix
%
%   with  A ~ U S V'
%
%-----------------------------------------------------------------------
%
%   D.C. Sorensen
%   10 Sep. 2013
%
%
%   Random search phase (Stage A)
%
     [m,n] = size(A);
     W = randn(n,2*k);


%
%    Get a basis for a 2k dim subspace of Range(A)
%
     Y = A*W;

%
%    Do q  power iterations with A A'
%
     for j = 1:q;

         T = A'*Y;
         Y = A*T;
         if (orth), [Y,G] = qr(Y,0); end  % Orthogonalize columns of Y

     end

%
%    Get an orthogonal basis U for the dominant subspace
%    of dimension 2k
%

     [U,R] = qr(Y,0);
%
%    Low Rank Factorization phase (Stage B)
%
%
%    Form matrix B with dominant singular value information
%
     B = U'*A;

     [U1,S,V] = svd(B);

     U = U*U1;

