function [X,R,H,Q]=FOM(A,b,x0);
% FOM full orthogonalization method
%   [X,R,H,Q]=FOM(A,b,x0) computes the decomposition A=QHQ', 
%   Q orthogonal and H upper Hessenberg, Q(:,1)=r/norm(r), 
%   using Arnoldi in order to solve the system Ax=b with the 
%   full orthogonalization method. X contains the iterates 
%   and R the residuals

n=length(A); X=x0;
r=b-A*x0; R=r; r0norm=norm(r); 
Q(:,1)=r/r0norm; 
for k=1:n
  v =A*Q(:,k);
  for j=1:k
    H(j,k)=Q(:,j)'*v; v=v-H(j,k)*Q(:,j);
  end
  e0=zeros(k,1); e0(1)=r0norm;                    % solve system
  y=H\e0; x= x0+Q*y;
  X=[X x];
  R=[R b-A*x];
  if k<n
    H(k+1,k)=norm(v); Q(:,k+1)=v/H(k+1,k);
  end
end