function PR = PageRank(G,alpha)
% G is an n-by-n link matrix in sparse format
% alpha is a positive parameter
% PR is permutation of 1:n with the property that webpage i
% has page rank PR(i) according to the Google model

[n,n] = size(G);
PR = zeros(n,1);

% Initializations
% Where are the zero and nonzero G-columns?
possums = zeros(1,n);
for j=1:n
possums(j) = sum(G(:,j));
end
ip = find(possums > 0);
iz = find(possums == 0);
% Scale the nonzero columns of G...
for k=ip
G(:,k) = (alpha/possums(k))*G(:,k);
end
% To minimize flops in the matrix-vector products, we precompute a
% pair of scaled vectors...
f = ones(n,1)/n;
% Initial guess will be f. We relax the termination criteria a bit
% noting that if x is the stationary vector and sum(x) = 1, then 1/n
% is the average spacing between the x(i).
opts = struct('disp',0,'v0',f,'tol',.001/n);
% Compute the stationary vector of the Google matrix..
[x,s,flag] = eigs(@(x)MyProd(x,G,f,alpha,ip,iz),n,1,'LM',opts);
% Page rank is based on the size of its components...
[z,idx] = sort(abs(x),'descend');
PR(idx) = 1:n;

function y = MyProd(x,G,f,alpha,ip,iz)
% y = (GoogleMatrix)*x
y = G(:,ip)*x(ip) + (sum(x(iz)) + sum(x(ip))*(1-alpha))*f;
% About 3n + 2*NZ flops where NZ = number of nonzeros in G