function [lambda, v] = rayleigh(A, vin, tol, kmax);
% RAYLEIGH  Rayleigh quotient iteration to find an eigenvalue
%           of a square matrix.  No guarantees.
% requires: mylu_check.m, myslash_check.m
% example: >> rayleigh([2 1 1; 0 2 2; 1 -1 4], [1 2 3]')

if nargin < 4, kmax = 20; end
if nargin < 3, tol = 1.0e-12; 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 - lambda*eye(m)) \ v
  [w, flag] = myslash_check(A - lambda*eye(m),v);
  % if LU decomp saw divide-by-zero then stop with current answer
  if flag==true, break, end  
  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