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;
