function [theta x iters]=powereig(A,maxiter)

% Power Method for estimating the dominant eigenvalue and corresponding
% eigenvector.  Only works if the dominant eigenvalue of A is distinct
% Thus, complex eigenvalues, which come in conjugate pairs, can not be
% estimated using this version of the power method (ie need a block
% version).
% 



n=size(A,1);
if n~=size(A,2)
    error('Matrix must be square')
end


if ~exist('maxiter','var')
    maxiter=10*n;
end

iters.x=zeros(n,maxiter);
iters.x(:,1)=randn(n,1);
iters.x(:,1)=randn(n,1)/norm(iters.x(:,1));

iters.theta=zeros(maxiter,1);
iters.normx=zeros(maxiter+1,1);
iters.normx(1)=1;

% Only if we do the ratios below
%[W eigA]=eig(A);
%eigA=diag(eigA)


for i=1:maxiter
    y=A*iters.x(:,i);
    
    iters.theta(i)=iters.x(:,i)'*y;   % Rayleigh Quotient (x'Ax)/(x'x)
    
    % Get theta this way if x is scaled by something other than its norm
    % iters.theta(i+1)=iters.x(:,i+1)'*y/(iters.x(:,i+1)'*iters.x(:,i+1)); 
    
    iters.x(:,i+1)=y/norm(y);
    
    % This is how Watkins p315 scales each estimated eigenvaector
    %iters.x(:,i+1)=y/norm(y,inf);
    
    
    
    % This ratio converges to (next to largest eig)/(largest eig) when
    % x is scaled by infinity norm (Watkins p.316)
    %if norm(iters.x(:,i)-W(:,end)/W(1,end))>0
    %  iters.ratio(i+1)=norm(iters.x(:,i+1)-W(:,end)/W(1,end))/norm(iters.x(:,i)-W(:,end)/W(1,end));
    %end
    
    % This ratio converges to 1 when x is scaled by Euclidean norm
    %if norm(iters.x(:,i)-W(:,end)/W(1,end))>0
    %  iters.ratio(i+1)=norm(iters.x(:,i+1)-W(:,end))/norm(iters.x(:,i)-W(:,end));
    %end
    iters.normx(i+1)=norm(iters.x(:,i+1)-iters.x(:,i));
end

theta=iters.theta(end);
x=iters.x(:,end);


figure(1)
clf
plot(abs(iters.theta))
title('\theta_i')

figure(2)
clf
semilogy(abs(diff(iters.theta)))
title('Check convergence: |\theta_{i+1}-\theta_i|')

figure(3)
clf
semilogy(iters.normx)
title('Check convergence: ||x_{i+1}-x_i||')


%plot(iters.ratio)
%title('Check convergence: ||x_{i+1} - w||/||x_i - w||')

% For comparisons
eigA=eig(A)
abs(eigA)