function hw7p1

rng(0)
m     = 25;
b     = 5;
A     = randn(m);
A     = A + A';
niter = 30;

%%% Compute the "exact" eigenvalues and order them by magnitude:
ee    = eig(A);
[~,ind] = sort(-abs(ee));
ee = ee(ind);

%%% Run the power iteration.
L_EST = zeros(niter,b);
B = randn(m,b);
for k = 1:niter
  B          = A*B;
  
  %%% Compute the evals of A restricted to the range of B:
  [Q,~]      = qr(B,0);
  eeloc      = eig(Q'*A*Q);
  [~,ind]    = sort(-abs(eeloc));
  L_EST(k,:) = eeloc(ind)';
end


kk  = 1:niter;
ERR = abs(L_EST - ones(niter,1)*ee(1:b)');
semilogy(kk',ERR)
legend('Convergence of lambda_1','Convergence of lambda_2','Convergence of lambda_3')

keyboard
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
