function [V,H,T,f] = RKA(E,A,k,v,s);
%
%   This function computes an Arnoldi-like 
%   factorization giving an ortho-normal basis
%   for a rational Krylov space (defined by A,E,s)
%   through and implementation of the RKA of Ruhe.
%
%   Input:  E -- an n by n matrix
%           A -- an n by n matrix
%           k -- a positive integer (k << n assumed)
%           v -- an n  vector (v .ne. 0 assumed)
%           s -- a vector holding k "shifts"
%
%   Output: V -- an n by k orthogonal matrix
%           H -- a k by k upper Hessenberg matrix
%           T -- a k by k upper triangular matrix
%           f -- an n vector
%
% 
%           with   AVH = EV(HS + T) - (A - s(k) E) f e_k' 
%
%           where S = diag(s)
%
%           assures V'V = I_k  and  V' f = 0 with  DGKS correction
%
%
%   D.C. Sorensen
%   17 Feb 09
%
    n = length(v);
    H = zeros(k);
    T = zeros(k);
    V = zeros(n,k);

    v = v/norm(v) ;

    w = (A - s(1)*E)\E*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; T(1,1) = 1;
        

    for j = 2:k,
       
        beta = norm(f);
	v = f/beta;

        H(j,j-1) = beta; 
        V(:,j)   = v;

	%   
	% The following choice of t  is arbitrary
	% and should be investigated further
	%
	t = rand(j,1);

        w = (A - s(j)*E)\E*V(:,1:j)*t;
        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; T(1:j,j) = t;

    end