function [ ret ] = fastpr( A , mode , alpha)
%FASTPR   Calculates simple Google PageRank using the adj. matrix
%
%  OUT = FASTPR( A , mode , alpha) where A is the adjacency matrix 
%    of the network.  mode  results in the output to be 
%   - 'p' : PageRanks (unsorted)
%   - 'po': PageRanks (sorted in decr. order)
%   - 'id': Docids (sorted in decr. order of importance) (default)
%
%    alpha  is optional, defaulted to 0.85
%
%  Last change:  10/2/2005 - Florian Knorn

if nargin < 1 || nargin > 3
    error('Please input adjacency matrix (and mode) (and alpha) !');
else
	N = size(A,1);
	if N ~= size(A,2)
		error('Please input SQUARE matrix !');
	end
	if nargin <= 2
		alpha = 0.85;
	end, if nargin <= 1
		mode = 'id';
	end
end

err = 42;							% set initial error to some value
eps = 10^(-floor(log10(N)+1)-3); % set error tolerance

P = double(A); % make sure entries of P are doubles (and not bools)

one_over_ns = ones(1,N) / N;	% n element vector with 1/n--entries
p = one_over_ns;					% initial PR distribution

% 1. manipulate matrix: adj -> raw google matrix
rowsums = sum(A,2);				% calc all the row-sums
zrows = find(rowsums==0);		% get indices of all the zero-rows

% replaces zeros by ones to prevent division by zero in next step
rowsums(zrows) = ones(length(zrows),1); 

P = spdiags(1./rowsums,0,N,N)*P; % devide each row by its rowsum
P(zrows,:) = ones(length(zrows),N)/N; % fix dangling nodes

% 2. perform the power method
P = alpha*P; z = (1-alpha)*one_over_ns; i = 1;
while err > eps 
	i = i+1; p_old = p;			% increase counter, store old p
	p = p_old*P + z;				% update p
	err = max( abs(p-p_old) );	% error
end

% return results

if strcmp(mode,'p') && ~strcmp(mode,'op') && ~strcmp(mode,'po')
	ret = p;
elseif strcmp(mode,'po') || strcmp(mode,'op')
	ret = sort(p, 'descend');
elseif strcmp(mode,'id')
	[p,i] = sort(p, 'descend');
	ret = i;
else
	error('unknown mode !');
end