function [x,rhovec] = myGMRES(A,b,k,tol)
% A naive implementation of GMRES
% to solve Ax = b


rhovec = [];

Q(:,1) = b/norm(b);

for j = 1:k
        
    v = A*Q(:,j);
    
    for l = 1:j
        H(l,j) = Q(:,l)'*v;
        v = v - H(l,j)*Q(:,l);
    end
    
    H(j+1,j) = norm(v);
    Q(:,j+1) = v/H(j+1,j);
    
    % After j steps, Q_{j+1} and H_{j+1} are formed
    % Solve the LSP
    Qb = Q'*b;
%    keyboard
    alpha = H\(Qb);
    
    % Form the residual and the approximate solution
    rho = norm(Qb - H*alpha);
    x = Q(:,1:j)*alpha;
    
%    keyboard
    
    rhovec = [rhovec log10(rho)];
    
    if (rho < tol)
        return;
    end
end

return;

