function [H,Q,v]=Arnoldi(A,b,k)
% ARNOLDI Arnoldi iteration
%   [H,Q,v]=Arnoldi(A,b,k) applies k<=n steps of the Arnoldi iteration to
%   the matrix A starting with the vector b. Computes Q orthogonal and
%   H upper Hessenberg such that AQ=QH+ve_k^T, with Q(:,1)=b/norm(b).

Q(:,1)=b/norm(b);
for j=1:k
  v=A*Q(:,j);
  for i=1:j
    H(i,j)=Q(:,i)'*v;
    v=v-H(i,j)*Q(:,i);
  end
  if j<k
    H(j+1,j)=norm(v);
    if H(j+1,j)<1e-15, disp('lucky breakdown'); break; end
    Q(:,j+1)=v/H(j+1,j);
  end
end
