function [v, lambda] = eig_power(A, x0, max_iter, delta)
% [v, lambda = eig_power(A, x0, max_inter=100, delta=1e-10)
%
% Find the largest eigenvalue and the corresponding eigenvector
% using the power iterations method.
%
% Parameters:
%    A : a dense square matrix (with real elements) [n x n]
%    x0: initial guess for the eigenvector. If not provided,
%        it will be taken as a random permutation of 1:n divided
%        by n.
%    max_iter: (default 100) maximum number of iterations
%    delta: (default 1e-10) if two successive iterations do not 
%        change the solution by more than delta, then stop
% Return:
%    lambda: eigenvalue
%    v: eigenvector

[m,n] = size(A);
if m ~= n
    error('The matrix A must be square!');
end
if nargin < 2
    x0 = randperm(n) ./ n;
end
if length(x0) ~= n
    error('Initial vector x0 must have the same size as the columns of A!');
end
if nargin < 3
    max_iter = 100;
end
if nargin < 4
    delta = 1e-10;
end

x0 = x0(:);                  % make it a column vector
lambda = 0;

for k = 1:max_iter
    lambda1 = lambda;
    x1 = A*x0/norm(x0,inf); 
    [xm,m] = max(abs(x1));
    lambda = x1(m);          % new eigenvalue: the largest in new x
    
    if (norm(x1-x0) < delta) && (abs(lambda1-lambda) < delta)
        break
    end
    
    x0 = x1;                 % new approximation is x0
end

if k == max_iter
    warning('Iterations exhausted before convergence!');
end

v = x0/norm(x0);             % ensure unit norm

return
