%%%%%%%%%%%%%% qr1.m
% QR Movie
% No Hessenberg reduction, no shifts

 figure(1), clf
 set(gcf,'doublebuffer','on')
 N = 20;
 V = randn(N);              % random eigenvector matrix
 A = V*diag([1:N])*inv(V);  % A has eigenvalues 1, 2, ..., 25, random eigenvectors
% A = randn(N);              % contrast with the previous A
 
 I = eye(N); j=0;
 x = kron([1:N+1]-.5,[1 1]); x = x(2:end-1);
 A = hess(A);
 [Q,R] = qr(A); A = R*Q; j=j+1;
 pcolor(x,x,kron(log10(abs(A)+1e-16),[1 1;1 1]));
 axis equal, axis([.5 N+.5 .5 N+.5])
 set(gca,'Ydir','Reverse') 
 caxis([-16,1]), colorbar, hold on
 shift = N;
 while (norm(diag(A,-1)) > 1e-14)
    if (shift>1) & (abs(A(shift,shift-1))<1e-20), shift=shift-1; end
    mu = A(shift,shift);
    X = A(1:shift,1:shift)-mu*eye(shift);
    [Q,R] = qr(X(1:shift,1:shift)); A(1:shift,1:shift) = R*Q+mu*eye(shift); j=j+1;
    pcolor(x,x,kron(log10(abs(A)+1e-16),[1 1;1 1]));
    axis equal, axis([.5 N+.5 .5 N+.5])
    set(gca,'Ydir','Reverse') 
    title(sprintf('Basic QR Algorithm: iteration %d', j),'fontsize', 16);
%    text(N+(N/20)*2.2,N+(N/20),'color indicates')
%    text(N+(N/20)*2.7,N+(N/20)*2.1,'log_{10} | a_{ij}^{(k)} |')
    text(N+(N/20)*1.8,N+N/15,'color indicates')
    text(N+(N/20)*2.0,N+(N/15)*2,'log_{10} | a_{ij}^{(k)} |')
%    drawnow; 
    pause(0.5);  % uncomment and adjust to slow the display down 
 end;