function [v,rayquo,resid] = powermethod(A,v,iter)
% basic iteration of the power method: see A&G p.222
for j=1:iter % should replace this by a while loop that checks both
           % the max number of iterations and the norm of the 
           % residual r = A*v - lambda*v
    v = A*v; % matriv-vector product
    v = v/norm(v);
    % Rayleigh quotient: if v is an eigenvector, rayquo(j) is an eigenvalue
    % The Rayleigh quotient is v'*A*v/(v'*v), but in our case, v'*v=1
    rayquo(j) = v'*A*v; % this is inefficient because we should not
                     % compute A*v twice inside the loop, only once
    resid(j) = norm(A*v - rayquo(j)*v); % again, should not compute A*v again!
end