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

% QUADPAGERANK  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
%               quadratic extrapolation applied every l iterations.
%
% EXAMPLE: [pi,time,numiter]=quadpagerank(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 = quadratic 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 
    % 'quadratic extrapolation'
      nextpi=alpha*pi*H + (alpha*(pi*a)+1-alpha)*v;
      nextnextpi=alpha*nextpi*H + (alpha*(nextpi*a)+1-alpha)*v;
      y=pi-prevpi;   nexty=nextpi-prevpi;   nextnexty=nextnextpi-prevpi;
      Y=[y' nexty'];
      gamma3=1;
      % do modified gram-schmidt QR instead of matlab's [Q,R]=qr(Y);
      [m, n] = size(Y);
      Q = zeros(m,n);
      R = zeros(n);
      for j=1:n
          R(j,j) = norm(Y(:,j));
          Q(:,j) = Y(:,j)/R(j,j);
          R(j,j+1:n) = Q(:,j)'*Y(:,j+1:n);
          Y(:,j+1:n) = Y(:,j+1:n) - Q(:,j)*R(j,j+1:n);
      end
      Qnextnexty=Q'*nextnexty';
      gamma2=-Qnextnexty(2)/R(2,2);
      gamma1=(-Qnextnexty(1)-gamma2*R(1,2))/R(1,1);
      gamma0=-(gamma1+gamma2+gamma3);
      beta0=gamma1+gamma2+gamma3;
      beta1=gamma2+gamma3;
      beta2=gamma3;
      nextnextpi=beta0*pi+beta1*nextpi+beta2*nextnextpi;
      nextnextpi=nextnextpi/sum(nextnextpi);
      pi=nextnextpi;
      %'end quadratic extrapolation'
  end
  pi=pi/sum(pi);
end
numiter=k;
time=toc;