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

ninf=norm(A,inf);
Q(:,1)=r/norm(r);
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 abs(H(j+1,j))+ninf==ninf
      disp('lucky breakdown occured!');
      H=H(1:end-1,:);
      break;
    end
    Q(:,j+1)=v/H(j+1,j);
  end
end