function [u,k] = generalized_minres(A,b,u0,epsilon)

n = length(u0);

r0 = A*u0-b;
rho0 = norm(r0);

v0 = r0/rho0;

V = zeros(n,n);
V(:,1) = v0;

R = zeros(n,n);

c = rho0 * [1;zeros(n-1,1)];

a = zeros(n,1);
b = zeros(n,1);

for k=1:n
    
    w = A*V(:,k);
    
    vHat = w;
    for l=1:k
        beta = V(:,l)'*w;
        R(l,k) = beta;
        vHat = vHat - beta*V(:,l);
    end
    beta = norm(vHat);
    
    if beta < 1e-8
        for j=1:k-1
            Gj = [ a(j) , b(j) ; -b(j) , a(j) ];
            R(j:j+1,k) = Gj * R(j:j+1,k);
        end    
        break;
    end
    
    R(k+1,k) = beta;
    V(:,k+1) = vHat/beta;
    
    for j=1:k-1
        Gj = [ a(j) , b(j) ; -b(j) , a(j) ];
        R(j:j+1,k) = Gj * R(j:j+1,k);
    end    
    
    denom = sqrt(R(k,k)^2+R(k+1,k)^2);
    a(k) = R(k,k)/denom;
    b(k) = R(k+1,k)/denom;
    Gk = [ a(k) , b(k) ; -b(k) , a(k) ];
    
    R(k:k+1,k) = Gk * R(k:k+1,k);
    
    c(k:k+1) = Gk*c(k:k+1);
    rho = abs(c(k+1));
        
    rho/rho0
    if rho/rho0 < epsilon
        break;
    end
    
end

alpha = zeros(k,1);
for j=k:-1:1
    alpha(j) = c(j)/R(j,j);
    c(1:j-1) = c(1:j-1) - R(1:j-1,j)*alpha(j);
end

u = u0 - V(:,1:k)*alpha;

