%  Demo_power.m
%
%  Demonstrate power method to compute maximal eigenpair of matrix A.

  %  Enter matrix A and compute dominant eigenvalue and corresponding 
  %  eigenvector. Also, get ratio 2 largest eigenvalues.

%A = [2 -1 0; -1 2 -1; 0 -1 2];  %  matrix A
A = [1 2; 3 4];  %  matrix A
[v,d] = eig(A);  %  eigenvalues and eigenvectors of A
d = diag(d);  %  Diagonal entries of d are eigenvalues of A.
[lambda,indx] = sort(d,'descend');  %  Sort eigenvalues in descending order.
lambda_max = d(indx(1));  %  dominant eigenvalue of A
v_max = v(:,indx(1));     %  corresponding eigenvector of A
ratio = abs( d(indx(2))/d(indx(1)) );  %  ratio of 2 largest eigenvalues

  [m,n] = size(A);  %  Get size of matrix.
  if m ~= n
    error(' *** A must be a square matrix ***');
  end

  max_iter = input(' Max. no. of iterations = ');

  %  Initialize.

  e_vec = [];
  theta_vec = [];
  q = randn(n,1);  %  Initial guess with iid Gaussian entries.
  u = q / norm(q)  %  Normalize initial guess.

  if n == 2  %  Plot estimate for maximal eigenvector.
    figure(1)  %  Put plot in figure 1.
    plot([0,u(1)], [0,u(2)], 'o-')
    xlabel('1st component')
    ylabel('2nd component')
    title('Normalized Eigenvector')
    axis('equal'), hold on  %  This "holds" plot in figure 1.
  end

  %  Iterate

  for k=1:max_iter
    q = A*u;
    beta = q'*u      %  Estimate for dominant eigenvalue.
    u = q / norm(q)  %  Estimate for corresponding eigenvector.
    fprintf(' Hit any key to continue.\n');
    pause

    if n == 2  %  Plot estimate for maximal eigenvector.
      plot([0,u(1)], [0,u(2)], 'o-'), axis('equal')
    end

    %  Compute error indicators.

    e_k = abs(beta - lambda_max);  %  eigenvalue error indicator.
    theta_k = acos( abs( u'*v_max ) );  %  eigenvector error indicator.

    %  Save error indicators in vectors e_vec and theta_vec.

    e_vec = [e_vec; e_k];
    theta_vec = [theta_vec; theta_k];

  end  %  end loop on k

  if n == 2, hold off; end  %  This "releases" plot in figure 1.

  %  Plot vectors that contain error indicators.

  figure(2)  %  Put plot in figure 2.
    semilogy(e_vec,'o-')
    xlabel('Iteration k')
    ylabel('|beta_k - \lambda_{max}|')
    title('Eigenvalue Error vs Power Method Iteration')
  figure(3)  %  Put plot in figure 3.
    semilogy(theta_vec,'o-')
    xlabel('Iteration k')
    ylabel('|angle between u_k and v_{max}|')
    title('Eigenvector Error vs Power Method Iteration')

  %  Compute and display convergence rate constants.

  slope1 = mean(diff(log(e_vec)));  %  Average slope in log e_k plot.
  c1 = exp(slope1);
  slope2 = mean(diff(log(theta_vec)));  %  Average slope in log theta_k plot.
  c2 = exp(slope2);

  fprintf(' Estimated eigenvalue convergence constant c1  = %6.4f\n', c1);
  fprintf(' Estimated eigenvector convergence constant c2 = %6.4f\n', c2);
  fprintf(' Abs value of ratio of 2 largest eigenvalues = %6.4f\n', ratio);