function [T,Z] = orth_iter(A, Z0, nsteps, displ)
%
% call: [T,Z] = orth_iter(A, Z0, nsteps, displ)
%
% Orthogonal iteration with A applied to block of vectors in Z0
% (also known as subspace iteration, simultaneous iteration, block power)
% Returns Z'AZ (which becomes triangular in limit under some assumptions)
% as well as Z.
%
Z = Z0;
for i = 1:nsteps
Y = A*Z;
[Z,U] = qr(Y,0);
if displ
Z'*A*Z
pause
end
end
T = Z'*A*Z;