%%
clear
clc
rng(8)

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

%%
% We use the Arnoldi iteration to generate orthonormal vectors. 
ndim = 60;
[Q, H] = arnoldi(A, b, ndim);

%% 
% The Arnoldi bases are used to solve the least squares problems defining
% the GMRES iterates. 
resid = zeros(ndim, 1);
nb = norm(b);
for m = 1:ndim  
    s = [nb; zeros(m, 1)];
    z = H(1:m+1, 1:m)\s;
    x = Q(:, 1:m)*z;
    resid(m) = norm(b - A*x);
 end

%%
% The approximations converge smoothly, practically all the way to machine epsilon.
semilogy(1:ndim, resid, '.-')
xlabel('m')
ylabel('|| b-Ax_m ||')
title('Residual for GMRES')
grid on