function [V,T] = blockpowermethod(A,V,iter)
% block power method beginning with V, which must have independent columns
% see A&G p.222, where it is known as subspace iteration
for j=1:iter % should control this by a while loop that checks both
           % the max number of iterations and the norm of some kind of residual 
    V = A*V; % matrix-vector product
    [V,R] = qr(V,0); % orthogonalize the columns of V: essential to make this work
    
    T = V'*A*V; % this is inefficient because we should not
                     % compute A*V twice inside the loop, only once
end