function [U,S,V] = randsvd(A, k, ell, ind)
    % 
    % Inputs 
    % A  : matrix of size n x p
    % k  : target rank, usually lower than the true rank
    % ell  : oversampling parameter; 
    % 
    % Outputs
    % U  : n x k (orthonormal columns)
    % S  : k x k (diagonal matrix)
    % V  : p x k (orthonormal columns)


    [n,p] = size(A);
    
    %Random sampling and sketching
    Omega = randn(p, k + ell);

    for j=1:ind
      ind1 = (j-1)*n+1;
      ind2 = j*n;
      Y(ind1:ind2,:) = A*Omega;
    end
    
    % Computing an orthonormal basis
    [Q,~] = qr(Y, 0);
    
    
    % 
    Bt = zeros(p,k+ell);
    for j=1:ind
      ind1 = (j-1)*n+1;
      ind2 = j*n;
      Bt = Bt + A'*Q(ind1:ind2,:);
    end
    B = Bt';
    [Ub,S,V] = svd(B, 'econ');
    
    U = Q*Ub;
    
    %Truncation
    U = U(:,1:k);
    S = S(1:k,1:k);
    V = V(:,1:k);
    
end