%
% Simple implementation of the FOM method for approximate
% solution of Ax=b,
% x0 = initial guess.
% m = number of vectors used to construct the Krylov subspace
%
% FOM is a projection method with 
% K=L= Krylov subspace K_m(A,r0)
% 
function x=fom(A,b,x0,m)
r0 = b - A*x0;
n = length(x0);
beta = norm(r0);
V = zeros(n,m);
V(:,1)=r0/beta;
H = zeros(m,m);
k=m;
for j=1:m
  w = A*V(:,j);
  for i=1:j
    H(i,j) = w'*V(:,i);
    w = w - H(i,j)*V(:,i);
  end
  tmp = norm(w);
  if tmp <= 1e-8
     k = j;
     break;
  end
  H(j+1,j) = tmp;
  V(:,j+1)=w/H(j+1,j);
end
z = beta*eye(k,1);
y = H(1:k,1:k)\z;
x = x0 + V(1:n,1:k)*y;
fprintf(1,'k=%d residual norm = %e\n',k,tmp*abs(y(k)) );

