function [V,H,f] = ArnoldiC(A,k,v);
%
%   Input:  A -- an n by n matrix
%           k -- a positive integer (k << n assumed)
%           v -- an n  vector (v .ne. 0 assumed)
%
%   Output: V -- an n by k orthogonal matrix
%           H -- a k by k upper Hessenberg matrix
%           f -- an n vector
%
% 
%           with   AV = VH + fe_k'
%
%           assures V'V = I_k with  DGKS correction
%
%
%   D.C. Sorensen
%   10 Feb 00
%
    n = length(v);
    H = zeros(k);
    V = zeros(n,k);

    v = v/norm(v);
    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;
        

    for j = 2:k,
       
        beta = norm(f);
%
%       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