% Matlab code qriter.m
% For "Applied Numerical Linear Algebra", Figures 4.3, 4.4
% Written by James Demmel, Oct 22, 1995
% Modified, Jun 5, 1997
%
% Illustrate convergence of unshifted QR iteration
% Used to produce figures 4.2 and 4.3
%
% Inputs:
% D = vector of eigenvalues
% m = number of iterations
%
% Outputs:
% plots of convergence of Orthogonal Iteration (or QR iteration)
% applied to a matrix with eigenvalues D and randomly chosen
% eigenvectors
%
[Ds,Is] = sort(-abs(D));
Ds = D(Is);
n = length(D);
Ds = reshape(Ds,n,1);
S = randn(n,n);
disp(['condition number(eigenvector matrix) = ',num2str(cond(S))])
A = S*diag(Ds)*inv(S);

Asave=A;
econv=0; clear econv;
rconv=0; clear rconv;
econv(1:n,1) = abs( diag(A) - Ds );
for j=1:n
rconv(j,1) = norm(A(j+1:n,1:j));
end
disp('Starting matrix='), Asave
disp('(pause - hit return to continue)'), pause
for i=1:m+1
[Q,R]=qr(A);
Alast = A;
A = R*Q;
disp(['Matrix ',int2str(i),' = ']),A
disp('(pause - hit return to continue)'), pause
econv(1:n,i) = abs( diag(A) + sort(-Ds) );
for j=1:n-1,
rconv(j,i) = norm(A(j+1:n,1:j));
end
end
figure(1), hold off, clf
figure(2), hold off, clf
sbplt = 0;
for i = 1:n,
if ( rem(sbplt,4) == 0 ) & ( i ~= 1 ),
disp('hit return to continue'), pause,
figure(1), hold off, clf
figure(2), hold off, clf
sbplt = 0;
end
sbplt = sbplt+1;
figure(1)
subplot(2,2,sbplt)
if i == 1,
rat = abs(Ds(2)/Ds(1));
elseif i == n,
rat = abs(Ds(n)/Ds(n-1));
else
rat = max(abs(Ds(i)/Ds(i-1)),abs(Ds(i+1)/Ds(i)));
end
dconv = (rat*ones(m+1,1)).^((0:m)');
semilogy(econv(i,:),'-b'), grid, hold on
semilogy(dconv,'--b'),
title(['Convergence of lambda(',int2str(i),')= ',num2str(Ds(i))])
xlabel('actual error (solid), estimate (dashed)')
if ( i < n )
figure(2)
subplot(2,2,sbplt)
rrat = abs(Ds(i+1)/Ds(i));
% rconv = (rrat*ones(m+1,1)).^((0:m)');
% semilogy(econv(i,:),'-b'), grid, hold on
% semilogy(dconv,'--b'),
dconv = (rrat*ones(m+1,1)).^((0:m)');
semilogy(rconv(i,:),'-b'), grid, hold on
semilogy(dconv,'--b'),
title(['Convergence of Schur form(',int2str(i+1),':',int2str(n),' , ', ...
'1:',int2str(i),')']),
xlabel('actual error (solid), estimate (dashed)')
end
end