function [V,H] = arnoldi(A,v,k)
%ARNOLDI  Generates Arnoldi factorisation of lenght k.
%   Given an n x n matrix A, an n x 1 vector v, and
%   a specified lenght k, with 1 <= k <= n-1,
%   [V,H] = arnoldi(A,v,k) generates an orthonormal
%   basis V_k+1 of the Krylov subspace K_k+1(A,v), 
%   together with the upper Hessenberg representation
%   H_k+1,k of A on this basis.
%
%   The following will hold:
%   A * V_k = V_k+1 * H_k+1,k
%
%   V is n x k+1 with orthonormal columns.
%   H is k+1 x k upper Hessenberg.
%   

% Use the start vector v as first column
V = v/norm(v);
H = [];

% Loop for successively larger subspaces
for i = 1:k
    w = A*V(:,i);               % Apply A to last column
    h = V'*w;                   % Project on previous columns
    w = w - V*h;                % Orthogonalize
    n = norm(w);                % Find norm
    w = w/n;                    % Normalize
    V = [V,w];                  % Add as new last column
    H = [H,h;zeros(1,i-1),n];   % Store coefficients
end