function [pi,time,numiter]=aitkenpagerank(pi0,H,v,n,alpha,epsilon,l);

% AITKENPAGERANK  computes the PageRank vector for an n-by-n Markov  
%                 matrix H with starting vector pi0 (a row vector), 
%                 scaling parameter alpha (scalar), and teleportation
%                 vector v (a row vector).  Uses power method with
%                 Aitken extrapolation applied every l iterations.
%
% EXAMPLE: [pi,time,numiter]=aitkenpagerank(pi0,H,v,900,.9,1e-8,10);
%
% INPUT: 	pi0 = starting vector at iteration 0 (a row vector)
%         H = row-normalized hyperlink matrix (n-by-n sparse matrix)
%         v = teleportation vector (1-by-n row vector)
%         n = size of P matrix (scalar)
%         alpha = scaling parameter in PageRank model (scalar)
%         epsilon = convergence tolerance (scalar, e.g. 1e-8)
%         l = Aitken extrapolation applied every l iterations (scalar)
%
% OUTPUT:   pi = PageRank vector 
%           time = time required to compute PageRank vector
%           numiter = number of iterations until convergence
%
% The starting vector is usually set to the uniform vector, 
% pi0=1/n*ones(1,n).
% NOTE: Matlab stores sparse matrices by columns, so it is faster  
%       to do some operations on H', the transpose of H.


% get "a" vector, where a(i)=1, if row i is dangling node 
%   and 0, o.w.

rowsumvector=ones(1,n)*H';
nonzerorows=find(rowsumvector);
zerorows=setdiff(1:n,nonzerorows); l=length(zerorows);
a=sparse(zerorows,ones(l,1),ones(l,1),n,1);


k=0;
residual=1;
pi=pi0;
tic;

while (residual >= epsilon) 
  prevpi=pi;
  k=k+1;
  pi=alpha*pi*H + (alpha*(pi*a)+1-alpha)*v;
  residual=norm(pi-prevpi,1);
  if (mod(k,l))==0 
    % 'Aitken extrapolation'
      nextpi=alpha*pi*H + (alpha*(pi*a)+1-alpha)*v;
      g=(pi-prevpi).^2;
      h=nextpi-2*pi+prevpi;
      nextpi=prevpi-(g./h);
      if (any(nextpi==-Inf)==1)
          pi=pi;
        else
            pi=nextpi;
      end
      %'end Aitken extrapolation'
  end
end
numiter=k;
time=toc;