% Matlab code powerplot.m
% For "Applied Numerical Linear Algebra", Algorithm 4.1
% Written by James Demmel, Oct 23, 1995
%                Modified, Jun  5, 1997
% 
% Plot convergence of the power method, Algorithm 4.1
% function [x,vv,verror,eerror] = powerplot(a,v,m)
%
% Inputs:
%  a = matrix
%  v = starting vector
%  m = number of steps
%
% Outputs:
%  Plot of actual and predicted error versus iteration
%  x = result of 30+m steps of power method
%  vv = matrix of all intermediate results of power method
%  verror = difference between x and each column of vv
%  eerror = difference between x'*a*x and each vv'*a*vv
%
function [x,vv,verror,eerror] = powerplot(a,v,m)
vv = v/norm(v);
for i = 1:m
  vv(:,i+1)=a*vv(:,i);
  vv(:,i+1)= vv(:,i+1)/norm(vv(:,i+1));
end
x=vv(:,m+1);
for i = 1:2*m+30,
  x = a*x;
  x = x/norm(x);
end
vv=vv*diag(sign(vv(1,:)));
x=x*sign(x(1));

% compute error = angle between vectors
% (note that it does not converge to macheps!)
% verror = acos(x'*vv);

% compute error = sine of angle between vectors
% (note that it does not converge to macheps!)
% cosine = x'*vv;
% verror = sqrt((1- cosine).*(1+cosine)) ;

% compute error = norm of difference of vectors
verror = (max(abs(vv-x*ones(1,m+1))))';

figure(1)
hold off, clf
% plot error in approximate eigenvectors
subplot(211)
semilogy(verror,'-b')
title('Eigenvector error, blue = actual, red = predicted')
xlabel('Iteration number')
grid
hold on
e=sort(abs(eig(a)));
[n,n]=size(a);
rat=e(n-1)/e(n);
rat=(rat*ones(1,m+1)).^(0:m);
semilogy(rat,'-r')
%
% plot error in approximate eigenvalues
subplot(212)
eerror = abs( ( diag(vv'*a*vv) - x'*a*x ) / (x'*a*x) );
semilogy(eerror,'-b')
title('Eigenvalue error, blue = actual, red = predicted')
xlabel('Iteration number')
grid
hold on
semilogy(rat,'-r')