% krylov.m -- Generalized Minimum Residual Method
% Written Fall 2005 by Eric Olson for Math/CS 466/666
%
% x=krylov(A,b,mmax,epsilon)
%
% Perform mmax iterations of GMRES method to solve A.  For more
% information see Algorithm 6.9 in Iterative Methods for Sparse 
% Linear Systems by Yousef Saad.

function x=krylov(A,b,mmax,tol)
    P=[]; b1=norm(b); e1=1;
    Q(:,1)=b/b1;
    for m=1:mmax
        t=A*Q(:,m);
        for j=1:m
           t=t-(Q(:,j)'*t)*Q(:,j);
        end
        Q(:,m+1)=t/norm(t);
        e1(m+1,1)=0.0;
        P(m+1,:)=0.0;
        t=Q'*A*Q(:,m);
        for j=1:m-1
            U(j,m)=P(:,j)'*t;
            t=t-U(j,m)*P(:,j);
        end
        U(m,m)=norm(t);
        P(:,m)=t/U(m,m);
        rnorm=norm((P*P(1,:)'-e1)*b1);
        disp(sprintf('norm(m=%d)=%g',m,rnorm));
        if rnorm<=tol
            break;
        end
    end
    z=(U\P(1,:)')*b1;
    x=Q(:,1:m)*z;
end
