function [Q,H] = arnoldi(A,u,m)
% ARNOLDI   Arnoldi iteration for Krylov subspaces.
% Input:
%   A    square matrix (n by n)
%   u    initial vector
%   m    number of iterations
% Output: 
%   Q    orthonormal basis of Krylov space (n by m+1)
%   H    upper Hessenberg matrix, A*Q(:,1:m)=Q*H (m+1 by m)

n = length(A);
Q = zeros(n,m+1);  
H = zeros(m+1,m);
Q(:,1) = u/norm(u);
for j = 1:m
  % Find the new direction that extends the Krylov subspace.
  v = A*Q(:,j);
  % Remove the projections onto the previous vectors.
  for i = 1:j
    H(i,j) = Q(:,i)'*v;
    v = v - H(i,j)*Q(:,i);
  end
  % Normalize and store the new basis vector.
  H(j+1,j) = norm(v);
  Q(:,j+1) = v/H(j+1,j);
end
