function [x,rho,xf,rhof,ritz,hritz] = Gmres_Fom(A,b,k);
%
%   Usage     [x,rho] = Gmres_Fom(A,b,k);
%
%   Compares GMRES and FOM
%
%   Solves  Ax = b
%
%   Input:  A     -- an n by n matrix
%
%
%           k    -- a positive integer 
% 
%
%   Output: x    -- an n  vector... approx solution Ax = b
%
%           rho  -- norm(b - Ax)
%
%           ritz --  Ritz values = roots of FOM polynomial
%       
%           hritz -- Harmonic Ritz values = roots of GMres polynomial
%
%
%   D.C. Sorensen
%   24 Oct 2008
%  
    e1 = eye(k,1);
    ek = zeros(k,1); ek(k) = 1;
    theta = norm(b);
    v = b/theta;
    n = length(v);

    [V,H,f] = ArnoldiC(A,k,v);

    Rnorm = norm(A*V - V*H - f*ek')

    beta = norm(f);
    if (k < n),
       [Q,R] = qr([H ; beta*ek']);
       rho = abs(Q(1,k+1)*theta);
    else
       [Q,R] = qr(H);
       rho = 0;
    end
    q = Q(1,1:k)'; 
    y = R(1:k,1:k)\q;
    x = V*y*theta;

    yf = (H\e1);
    xf = V*yf*theta;
    rhof = beta*abs(yf(k));
    
    g = ((H')\ek)*(beta^2);
    hritz = eig( H + g*ek');
    ritz  = eig(H);