function [V,H,x,flg] = LongCG(A,b,tol,x)
n = length(b);
V = zeros(n,n);
H = zeros(n,n);
flg = 0;
if (nargin <4)
    x = zeros(n,1);
end
res = b - A * x;
TOL = norm(A)*1e-12;
resnorm = zeros(n,1);
fucval  = zeros(n,1);
for t = 1:n
    V(:,t) = res;
    for s = 1:t-1
        Av = A*V(:,s);
        H(t,s) = sum(V(:,t).*Av)/sum(V(:,s).*Av);
        V(:,t) = V(:,t)-H(t,s)*V(:,s);
    end
    alpha = sum(res.*V(:,t))/sum(V(:,t).*(A*V(:,t)));
    x     = x + alpha * V(:,t);
    res   = b - A * x;
    resnorm(t) = norm(res);
    funcval(t) = sum(x.*(A*x))-2*sum(b.*x);
    if (norm(res)<=tol)
        V = V(:,1:t);
        H = H(1:t,1:t);
        H(abs(H)<TOL) = 0;
        H2 = H;
        for s = 1:t-1
            H2(s+1,s) = 0;
        end
        f1 = ['Conjugate Gradient Method Converged'];
        f2 = ['Matrix Dimension = ', num2str(n)];
        f3 = ['# of iterations = ', num2str(t), '; tolerance  = ', num2str(tol)];
        f4 = ['# of non-zeros in H below first subdiagonal = ', num2str(nnz(H2))];
        f  = strvcat(f1,f2,f3,f4);
        figure(1);
        clf;
        spy(H);
        hold on;
        title(f);
        hold off;
        g = ['Residual Norms'];
        figure(2);
        clf;
        title(g); 
        hold on;
        semilogy(resnorm,'b.-');
        hold off;
        h = ['Function Values'];
        figure(3);
        clf;
        title(h); 
        hold on;
        plot(funcval,'b.-');
        hold off;
        return;
    end
end
flg = 1;
H2 = H;
for s = 1:t-1
    H2(s+1,s) = 0;
end
f1 = ['Conjugate Gradient Method Converged'];
f2 = ['Matrix Dimension = ', num2str(n)];
f3 = ['# of iterations = ', num2str(t), '; tolerance  = ', num2str(tol)];
f4 = ['# of non-zeros in H below first subdiagonal = ', num2str(nnz(H2))];
f  = strvcat(f1,f2,f3,f4);
spy(H);
hold on;
title(f);
hold off;
g = ['Residual Norms'];
figure(2);
clf;
title(g); 
hold on;
semilogy(resnorm,'b.-');
hold off;
h = ['Function Values'];
figure(3);
clf;
title(h); 
hold on;
plot(funcval,'b.-');
hold off;