function [Q,H]=arnoldi(A,b,m)
% [Q,H]=arnoldi(A,b,m)
% A simple implementation of the Arnoldi method.
% The algorithm will return an Arnoldi "factorization":
%   Q*H(1:m+1,1:m)-A*Q(:,1:m)=0
% where Q is an orthogonal basis of the Krylov subspace
% and H a Hessenberg matrix.
%
% Example:
%  A=randn(100); b=randn(100,1);
%  m=10;
%  [Q,H]=arnoldi(A,b,m);
%  should_be_zero1=norm(Q*H-A*Q(:,1:m))
%  should_be_zero2=norm(Q'*Q-eye(m+1))
% 
% Exercise material for the course SF2524 (KTH - Elias Jarlebring)
    
    n=length(A);;
    Q=zeros(n,m+1);
    Q(:,1)=b/norm(b);
    for k=1:m
        w=A*Q(:,k); % Matrix-vector product with last element
        
        % Orthogonalize w against columns of Q (single GS)
        h=Q(:,1:k)'*w;  
        worth=w-Q(:,1:k)*h;        
        beta=norm(worth);  
        
        Q(:,k+1)=worth/beta; % normalize

        % Put Gram-Schmidt coefficients into H-matrix
        H(1:(k+1),k)=[h;beta]; 
    end
end