function x=gmresJD(A,Q,kappa,r,k)
%
% k-step GMRES to solve the correction equation for JD-method:
%
%   (I-Q*Q') (A-kappa I) (I-Q*Q') x = -r
%
%  Input
%            A      (n-by-n)
%            Q      (n-by-p) vector
%            kappa  scalar shift
%            r      (n-by-1) vector
%            k      (integer) number of GMRES steps
%
%  Output
%
%            x      (n-by-1) vector
%
%                                            B/L 7/26/2009
%
n=size(A,1);
V = zeros(n,k+1); H = zeros(k+1,k);

reorth = 1;

beta=norm(r); V(:,1)=-r/beta;
for j=1:k,
    % (I-Q*Q') (A-kappa I) (I-Q*Q') * V(:,j)
    vh = V(:,j);  
    vh=vh-Q*(Q'*vh);
    vh=A*vh-kappa*vh;
    vh=vh-Q*(Q'*vh); 
%   >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%   by MGS
    for i=1:j,
        H(i,j)=( V(:,i) )' * vh;
        vh = vh - H(i,j)*V(:,i);
    end
    if reorth == 1,
       for i=1:j,
           tmp=( V(:,i) )' * vh;
           vh = vh - tmp*V(:,i);
           H(i,j)=H(i,j)+tmp;
       end
    end
%   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    H(j+1,j)=norm(vh);        
    V(:,j+1)=vh/H(j+1,j);
end
%
% Least square  min || H w - beta *e1 ||

e1=zeros(k+1,1); e1(1)=beta;
w=H\e1; x=V(:,1:k)*w;
 
