function [ lambda, x, k ] = power_method( A, x0, tol, kmax );
% POWER_METHOD power method for the computation of the largest eigenvalue
% (in modulus) of the matrix A (\lambda_1). We assume that A is square,
% |\lambda_1| > |lambda_i| for i=2,...n, and \lambda_1 non zero
% Stopping criterion based on the relative difference of successive 
% iterates of the eigenvalue.
%  [ lambda, x ] = power_method( A, x0, tol, kmax )
%  Inputs: A      = matrix (n x n)
%          x0     = initial vector (n x 1)
%          tol    = tolerance for the stopping criterion
%          kmax   = maximum number of iterations
%  Output: lambda = computed (largest) eigenvalue
%          x      = computed eigenvector correspoding to lambda
%          k      = number of iterations
%

x = x0;
y = x / norm( x );
lambda = y' * A * y;
err = tol + 1;
k = 0;
lambda_old = lambda;

while ( err >= tol && k <= kmax )
   
    x = A * y;
    y = x / norm( x );
    lambda = y' * A * y;
    
    err = abs( lambda - lambda_old ) / abs( lambda );
    
    lambda_old = lambda;
    k = k + 1;
    
end

return