function [res,x] = fom (A, b, x0, m, tol) 
% function [res,x] = fom (A, b, x0, m, tol) 
% solves Ax = b by m steps of FOM -- bare-bone version
%%-----------------------------------------------
tolmac = 1.e-16;
n = size(A,1);
r = b - A*x0;
ro = norm(r,2);
rs(1) = ro ;
res(1) = ro ;
V(1:n,1) = r/ro; 
for i=1:m 
    i1 = i + 1  ;
    V(1:n,i1) = A*V(1:n,i)  ;
%------     modified GS  ;
    for j=1:i
        t = V(1:n,j)'*V(1:n,i1)  ;
        H(j,i) = t  ;
	V(1:n,i1) = V(1:n,i1) - t*V(1:n,j)  ;
    end
%% apply previous eliminations
    for k=2:i
        H(k,i) = H(k,i) - lll(k) * H(k-1,i) ;
    end
    if (i > 1) 
       rs(i) = - lll(i) *rs(i-1);
    end
%% test for convergence
    t = norm(V(1:n,i1),2)  ;
    res(i1) = abs(rs(i)*t/H(i,i));
    if (res(i) < tol*ro  ||  i==m) 
       y = triu(H) \ rs';
       x = x0 + V(:,1:i)*y ;
       break;
    end
    H(i1,i) = t  ;
    lll(i1) = t / H(i,i);
    if (t  ~= 0.0)  
        t = 1.0 / t  ;
	V(1:n,i1) = V(1:n,i1)*t  ;
    end
end