   function [v,res] = powerit(A, nit, v) 
%% [v,res] = powerit(A, nit, v) 
%% does nit steps of the power iteration 
%% starting with (input) v
%% 
if (nargin < 3) 
     v = randn(size(A,1),1);
end 
%%
 for i=1:nit
     v = A*v;
     [t,k] = max(abs(v));
     alpha = v(k);
     if (alpha == 0) 
        error('ERR. division by 0 -- STOP') 
     else 
       v = v / alpha;
       res(i) = norm(A*v - alpha*v) ;
    fprintf(1,...
       ' iter %3d -> res=%12.4e alph=%12.4e\n',...
  i,res(i),alpha);
     end
 end

