function [lambda, y, its, norm_r] = potenzeOP (OP, n, x, tol, max_its)
%
% [lambda, y, its, norm_r] = potenzeOP (OP, n, x, tol, max_its)
%
% Calcola l'autovalore di modulo massimo dell'operatore OP. OP puo`
% essere, per esempio,:
%
% * OP = @(x) A*x;
% * [L, U, P] = lu (A); OP = @(x) U \ (L \ (P * x));
% * [L, U, P] = lu (A - mu * eye (4)); OP = @(x) U \ (L \ (P * x));
%
% n e` la dimensione di x.
if ((nargin <= 2) || isempty (x))
  x = randn (n, 1);
end
if ((nargin <= 3) || isempty (tol))
  tol = [1e-6, 1e-6];
end
if ((nargin <= 4) || isempty (max_its))
  max_its = 100;
end
if (isscalar (tol))
  tol = [tol, tol];
end
y = x / norm (x);
x = OP (y);
its = 1;
lambda = y' * x;
norm_r = norm (x - lambda * y);
while ((norm_r > tol(1) * abs (lambda) + tol(2)) && (its < max_its))
  y =  x / norm (x);
  x = OP (y);
  lambda = y' * x;
  its = its + 1;
  norm_r = norm (x - lambda * y);
end
if (norm_r > tol(1) * abs (lambda) + tol(2))
  warning('Impossibile raggiungere la tolleranza richiesta')
end
%!test
%! A = [-2, 1 ,0, 0;1, -2, 1, 0;0, 1, -2 ,1;0, 0, 0, 1];
%! OP = @(x) A * x;
%! lambda = potenzeOP (OP, 4);
%! assert (abs (lambda), max (abs (eig (A))), 1e-5)
%!test
%! A = [-2, 1 ,0, 0;1, -2, 1, 0;0, 1, -2 ,1;0, 0, 0, 1];
%! [L, U, P] = lu (A);
%! OP = @(x) U \ (L \ (P * x));
%! lambda = potenzeOP (OP, 4, [], [], 200);
%! assert (abs (1 / lambda), min (abs (eig (A))), 1e-5)
%!test
%! A = [-2, 1 ,0, 0;1, -2, 1, 0;0, 1, -2 ,1;0, 0, 0, 1];
%! mu = -0.5;
%! [L, U, P] = lu (A - mu * eye (4));
%! OP = @(x) U \ (L \ (P * x));
%! lambda = potenzeOP (OP, 4, [], [], 200);
%! ev = eig (A);
%! [~, idx] = min (abs (ev - mu));
%! assert (abs (mu + 1 / lambda), abs (ev(idx)), 1e-5)
