% Steve Penny
% AMSC 600
% 10/18/07
% HW 2_2c

west0479 = load('west0479.mat');
C = west0479.west0479;

n = length(C);
k = 100;

%%%%%%%%%%%%%%
% Initialize %
%%%%%%%%%%%%%%
Z = zeros(n,k);
W = zeros(n,k);
B = zeros(n,k);
Z(:,1) = ones(length(Z),1)/norm(ones(length(Z),1),2);


% iterate algorithm from 2b 
% to generate sigular values of C
for j=1:k
    
    % Part 1
    if j==1
        W(:,j) = C*Z(:,j);
    else
        W(:,j) = C*Z(:,j) - W(:,j-1)*B(j-1,j);
    end
    B(j,j) = norm(W(:,j),2);
    W(:,j) = W(:,j)/B(j,j);
    
    % Part 2
    Z(:,j+1) = C'*W(:,j) - Z(:,j)*B(j,j);
    B(j,j+1) = norm(Z(:,j+1),2);
    Z(:,j+1) = Z(:,j+1)/B(j,j+1);
    
end

Sig_C_approx = svd(B);

Sig_C = svd(full(C));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display Results:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Print some of the singular values
Sig_C(1:10)
Sig_C(10:20)
Sig_C(20:30)
Sig_C(end-10:end)

Sig_C_approx(1:10)
Sig_C_approx(10:20)
Sig_C_approx(20:30)
Sig_C_approx(end-10:end)

% plot the singular values
figure(1)
plot(Sig_C,Sig_C,'r.')
hold on
plot(Sig_C_approx,Sig_C_approx,'bo','MarkerSize',10.0)
hold off
ylabel('singular values')
xlabel('singular values')
set(gcf,'color',[1 1 1])
set(gca,'box','off')

% plot the upper singular values
figure(2)
upper_cut = 3*10^5;
plot(Sig_C(find(Sig_C >= upper_cut)),Sig_C(find(Sig_C >= upper_cut)),'r.')
hold on
plot(Sig_C_approx(find(Sig_C_approx >= upper_cut)),Sig_C_approx(find(Sig_C_approx >= upper_cut)),'bo','MarkerSize',10.0)
hold off
ylabel('singular values')
xlabel('singular values')
set(gcf,'color',[1 1 1])
set(gca,'box','off')

% plot the lower singular values
figure(3)
lower_cut = 0.5*10^5;
plot(Sig_C(find(Sig_C <= lower_cut)),Sig_C(find(Sig_C <= lower_cut)),'r.')
hold on
plot(Sig_C_approx(find(Sig_C_approx <= lower_cut)),Sig_C_approx(find(Sig_C_approx <= lower_cut)),'bo','MarkerSize',10.0)
hold off
ylabel('singular values')
xlabel('singular values')
set(gcf,'color',[1 1 1])
set(gca,'box','off')

return