function [lambda, v] = poweriter(A, vin, tol, kmax);
% POWERITER  Power iteration to find the largest eigenvalue
%            of a square matrix.  No guarantees.
% example: >> poweriter([2 1 1; 0 2 2; 1 -1 4], [1 2 3]')

if nargin < 4, kmax = 30; end
if nargin < 3, tol = 1.0e-4; end
if diff(size(A)) ~= 0, error('only square matrices have eigenvalues'), end
m = size(A,1);
if size(vin,1) ~= m, error('vin must be m x 1 if A is m x m'), end

v = vin / norm(vin);
lambda = v' * A * v;
if lambda == 0.0, return, end

for k = 1:kmax
  w = A * v;
  if norm(w)==0, lambda = 0; break, end % if divide by zero, stop with current answer
  v = w / norm(w);
  lamnew = v' * A * v;
  if abs(lambda - lamnew) < tol, break, end
  lambda = lamnew
end

if k>=kmax, error('max number of iterations exceeded'), end