% Before the call to eigs...
f = ones(n,1)/n;
g = (1-alfa)*f;
s = sum(G)/alfa;
zdx = find(s==0);   % The indices of pages with no outlinks
pdx = find(s>0);    % The indices of pages with at least one outlink;
% Scale the nonzero columns of G
for k=1:length(pdx)
    G(:,pdx(k)) = G(:,pdx(k))/s(pdx(k));
end




function y = MyProd(x,zdx,pdx,f,g,G,)
% G is the link matrix
% Carries out the matrix vector product y = Tx where T is the nxn defined by
%              e/n                  if sum(G(:,j)) = 0
%   T(:,j) = 
%              alfa*G(:,j)/sum(G(:,j)) + (1-alfa)*e/n
%
% where n = length(x); 
% Compute y as a linear combination of T's columns...
% Contributions from the no-outlink pages...
y = sum(x(zdx))*f;
% Now the contributions from the pages that have outlinks...
for j=1:length(pdx)
    k = pdx(j);
    col = g + G(:,pdx(j)); 
    y = y + x(k)*col;
end
