function [error,H] = fastsvdnew(A,s,k);
%
% [ROWS_PICKED,H] = fastsvd(A,s,k)
% Given an m-by-n matrix A, FASTSVD computes an approximation
% to the top k right singular vectors of A, by sampling s rows 
% of A. It also returns an m-vector ROWS_PICKED that contains
% 1 if a certain row was picked and 0 otherwise.  
%

% 
% Some initializations. The ER vector contains the norms of each 
% row of A. The PROB vector is built with the following formula:
% $PROB(i) = \sum_{j=1}^i ER(j)$.
%
[m,n]  = size(A);
NORMAFRO = norm(A,'fro')^2;
co = s/NORMAFRO;
rand('state',sum(100*clock));

for i=1:m,
 ER(i) = norm(A(i,:))^2;
end

PROB(1) = ER(1);
for i=2:m,
 PROB(i) = PROB(i-1) + ER(i);
end

%
% Pick a random subset of s rows of A to form S. We pick according
% to the probabilities described in our paper. The following proce-
% dure ensures it. We also normalize every row of A that is to 
% be included in S. This is also described in our paper.
%

% Initialize S to be an s-by-n all zeros matrix.

S = sparse(s,n);

for i=1:s,
  %
  % pick random number ROW from 1 up to NORMAFRO.
  %
  ROW = ceil (rand * NORMAFRO);    
  %
  % perform binary search in the PROB vector to locate the number 
  % ROW and pick the corresponding ROW of A to be included  in S.
  % 
  left = 1;                        
  right = m;                       
  while abs(left-right) > 1        
   if ROW > PROB(ceil((left+right)/2))  
      left = ceil((left+right)/2);
   else  right = ceil((left+right)/2); 
   end
  end
  %
  % We normalize the row of A to include it in S 
  %
  S(i,:) = A(left,:)/(sqrt(co*ER(left))); 
end

%
% compute the svd of S*S' (S' is the transpose of S). The 
% matrix G is a diagonal matrix containing the singular 
% values of S*S'. We need their square roots for a sub-
% sequent normalization. 
%
[P,G,Q]=svd(full(S*S'));
d = diag(G);
d = sqrt(d);

%
% divide the columns of P (right singular vectors - don't 
% forget that P = Q because SS' is symmetric) by the cor-
% responding (square root of) singular value. 
%
P = P(:,1:k);
i=1;
while i<=k 
 if abs(d(i)) >= 1e-2
  P(:,i) = P(:,i)/d(i);
 else
  P(:,i) = zeros(s,1);
 end
  i = i + 1;
end;

%
% form matrix H according to the paper
%
H = S'*P;

%
% Now we have computed the approximation that we need. 
% It will be the following matrix product: D = A*H*H'
% Here we compute the percentage of the Frobenius norm
% of A-AHH' (squared) over NORMAFRO. This is a measure 
% for the error (see paper).
% 

t = norm(A*H,'fro')^2;
error = 1 - t/NORMAFRO;