function hw7p1_extra

rng(0)
m     = 50;
b     = 5;
A     = fast_decay(m);
niter = 10;
ee    = eig(A);
[~,ind] = sort(-abs(ee));
ee = ee(ind);

%%% Run basic power iteration.
B_orig = randn(m,b);
B = B_orig;
for k = 1:niter
  B = A*B;
end

%%% Compute the eigenvalues resulting from the basic iteration.
[Q,~] = qr(B,0);
ee_basic = eig(Q'*A*Q);
[~,ind] = sort(-abs(ee_basic));
ee_basic = ee_basic(ind);

%%% Run stabilized power iteration from the same starting matrix.
Q = B_orig;
for k = 1:niter
  [Q,~] = qr(A*Q,0);
end

%%% Compute the eigenvalues resulting from the stabilized version.
ee_stab = eig(Q'*A*Q);
[~,ind] = sort(-abs(ee_stab));
ee_stab = ee_stab(ind);

format long
fprintf(1,'The exact eigenvalues:\n')
disp(ee(1:b)')
fprintf(1,'Approximations from basic method:\n')
disp(ee_basic(1:b)')
fprintf(1,'Approximations from stabilized method:\n')
disp(ee_stab(1:b)')
format short e
fprintf(1,'Errors from basic method:\n')
disp(ee(1:b)' - ee_basic(1:b)')
fprintf(1,'Errors from stabilized method:\n')
disp(ee(1:b)' - ee_stab(1:b)')
format

keyboard
return


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

function A = fast_decay(m)

xxs = rand(2,m);
xxt = [1.8+rand(1,m);rand(1,m)];
DD1 = xxs(1,:)'*ones(1,m) - ones(m,1)*xxt(1,:);
DD2 = xxs(2,:)'*ones(1,m) - ones(m,1)*xxt(2,:);
A   = log(DD1.*DD1 + DD2.*DD2);
A   = A + A';
A   = (1/norm(A))*A;

return

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