%%
clear
clc
rng(5)

%%
% First we define a triangular matrix with known eigenvalues and a random
% vector b.
np = 100;
lambda = 10 + (1:np);
A = diag(lambda) + triu(rand(np), 1); 
b = rand(np, 1);

%%
% Next we build up the first 30 Krylov matrices iteratively, using
% renormalization after each matrix-vector multiplication. 
ndim = 30;
Km = zeros(np, ndim);
Km(:, 1)  = b;
for m = 2:ndim      
    v = A*Km(:, m-1);
    Km(:, m) = v/norm(v);
end

%%
% Now we solve a least squares problem for Krylov matrices of increasing
% dimension.
warning off
resid = zeros(ndim, 1);
for m = 1:ndim  
    z = (A*Km(:, 1:m))\b;
    x = Km(:, 1:m)*z;
    resid(m) = norm(b - A*x);
end

%%
% The linear system approximations show smooth linear convergence at first,
% but the convergence stagnates after only a few digits have been found.
figure(1)
semilogy(resid, '.-')
xlabel('m')

ylabel('|| b - Ax_m ||')
set(gca, 'ytick', 10.^(-6:2:0))
title('Residual for linear systems')
grid on

%%
figure(2)
condition = zeros(ndim, 1);
for m = 1:ndim
    condition(m) = cond(A*Km(:, 1:m));
end
semilogy(1:ndim, condition, '.-')
xlabel('m')
ylabel('\kappa(m)')
title('Condition numbers of Krylov matrices')
grid on