function [V,H,f] = arnoldif(A,V,H,f,k,m)
%ARNOLDIF computes the Arnoldi decomposition of matrix A
%   Input:  A -- an n by n matrix
%           V -- an n x m orthogonal matrix
%           H -- an m x m upper Hessenberg matrix
%           f -- an n-dimensional column vector. on first call,
%                a random vector
%           k -- start of the Arnoldi steps
%           m -- dimension of the Krylov subspace, m << n.
%
%   Output: V -- an n by m orthogonal matrix
%           H -- an m x m upper Hessenberg matrix
%           f -- an n vector
%
%           AV = VH + f*e_m', where e_m is the standard basis
%           vector [0 0 ...0 1]'
%
%           assures V'V = I_m with  DGKS correction
%

%   This code was kindly contributed by
%   D.C. Sorensen
%   3 April 2013
%   Some reorganization was necessary in order to use it with eigsb

n = length(f);
if k == 1 
   beta = norm(f);
   v = f/beta;
   w = A*v;
   alpha = v'*w;

   f = w - v*alpha;
   c = v'*f;
   f = f - v*c;
   alpha = alpha + c;

   V(:,1) = v;
   H(1,1) = alpha;
   k = 2;
end
        

for j = k:m
       
   beta = norm(f);
   v = f/beta;
   w = A*v;
%
%  The following test asks if w = Av is numerically in range(V).
%  If the test is passed then beta is indistinguishable from 0
%  when compared to norm(w) and the component (I - VV')w is 
%  meaningless
%
   if ( beta < 10*eps*norm(w))
      %  Construct a random v and orthogonalize against V
      %  Since beta is indistinguishable from 0, v*beta will be of
      %  negligible size for any unit vector v
      %       
      v = randn(n,1);
      c = V'*v;
      v = v - V*c;
      c = V'*v;
      v = v - V*c;
      v = v/norm(v);
   else
      v = f/beta;
   end
        
   H(j,j-1) = beta; 
   V(:,j) = v;

   w = A*v;
   h = V(:,1:j)'*w;
   
   f = w - V(:,1:j)*h;
   c = V(:,1:j)'*f;
   f = f - V(:,1:j)*c;
   h = h + c;

   H(1:j,j) = h;
end 
