function [eigval,eigvec] = rayleigh_iter(A,x)
% Call : [lambda,v] = rayleigh(A) or
%	[lambda,v] = rayleigh(A,x) - x is the initial vector
%
% Task : Computes an eigenvalue and an associated
%	eigenvector of A using Rayleigh iteration.

[m,n] = size(A);

if m ~= n
	error('A must be square')
end

if nargin < 2
	x = randn(n,1);
	x = x/norm(x);
	rayleigh = x'*A*x;
end


epsln = 5*10^-6;
iter = 1;

fprintf('Rayleigh quotient qk*Aqk\n');

while iter == 1 | (abs(rayleigh-rayleighold) > epsln)
	xold = x;
	
	x = (A-rayleigh*eye(n))\x;
	x = x/norm(x);
	rayleighold = rayleigh;
	rayleigh = x'*A*x;
	
	
	fprintf('iter %d : %1.15f\n', iter, rayleigh);
	
	iter = iter+1;
	
end

eigvec = x;
eigval = rayleigh; 

return