clear ;

n = 5;
A = rand(n);
B = A'*A;
[Q R] = qr(B);
d = [7 4 1 9 2];
D = diag(d);
K = Q'*D*Q;

%option = 1; % Metodo das potencias
%option = 2; % Metodo das potencias inversas
option = 3 % Metodo das potencias inversas com deslocamento
q  = 3.4;

if(option == 1)
  M = K;
elseif(option == 2) 
  M = inv(K);
elseif(option == 3)
  M = inv(K- q*eye(n));
end

x = rand(n,1);
y = x / norm(x);
lam = y'*M*y;
for k=1:10000
  x = M*y;
  y = x / norm(x);
  lamnew = y'*M*y;
  if(abs(lamnew - lam) < 1e-14) 
   break ;
  end
  lam = lamnew;
end

if(option == 1)
 lam
elseif (option == 2)
 1/lam
elseif (option == 3)
 1/lam + q
end