function [x, err, iter, flag] = fastgmres( A, F1, G1, b, x, max_it, tol )
% [x, err, iter, flag] = fastgmres( A, F1, G1, b, x, max_it, tol )
% 
%  fastgmres.m solves the linear system Ax=b
%  using a fast version of the Generalized Minimal RESidual ( GMRES ) method
%  This fast version does not require restarting for symmetric + low rank matrices.  
%  We suppose that the matrix A-F1*G1' is hermitian. 
%  (c) Lothar Reichel and Bernhard Beckermann, August 21,2006
% 
%  input   A        real or complex nonsymmetric matrix
%          F1, G1   real or complex matrix such that A-F1*G1' is hermitian
%          b        real or complex right hand side vector
%          x        real or complex initial guess vector (optional)
%
%          max_it   INTEGER maximum number of iterations (optional)
%          tol      real  residual residual error tolerance (optional)
% 
%  output  x        real or complex solution vector
%          err      residual error norm (exact)
%          iter     INTEGER number of iterations performed
%          flag     INTEGER: 0 = solution found to tolerance
%                            1 = no convergence given max_it

% storage requirements (beside input A, F, G)
%    vectors of size of dimension of system:  
%            x, vzprime, v(j:j+1), z(j:j+1), Fbar(1:s), W(1:s)
%    vectors of size j: 
%            tridiagonal t (as sparse matrix=3 vectors)
%    vectors of size s:  
%            Fhat_jth_row, Ghat_jth_row, pstar 

% at beginning of iteration j, v=v_j, v_old=v_{j-1}, z=z_j, z_old=z_{j-1} 
% during iteration j: vzprime contains either vprime or zprime

n=size(A,1); 
s_half=size(F1,2);
% Notice that F=[F1,G1], G=[G1,-F1]
if nargin<4, error('We need at least 4 arguments');end;
if nargin<5, x=zeros(n,1);end;
if nargin<6, max_it=floor(n/3);end;
if nargin<7, tol=norm(b-A*x)*eps*1000;end;

flag=1;
j=1; 
%r=b-A*x;
v=b-A*x;
gamma=norm(v);
v=v/gamma;                    % v_1
t=spalloc(max_it+1,max_it+1,3*max_it);
c_angle=1; s_angle=0;         % initialization c_0, s_0
Fbar=zeros(n,2*s_half);       % initialisation f_{1:s,0} 
W=zeros(n,2*s_half);          % initialisation w_{1:s,0}
z=x/gamma;                    % z_1
% j=1
Fhat_jth_row=v'*[F1,G1];      % first row of Fhat
Ghat_jth_row=[Fhat_jth_row(s_half+1:2*s_half),-Fhat_jth_row(1:s_half)]; % first row of Ghat 
pstar=Fhat_jth_row;           % p1^*
Fbar=v*Fhat_jth_row;          % Fbar_1=f_{1:s,1} 
W=z*Fhat_jth_row;             % W_1 = w_{1:s,1} 

% new Arnoldi vector
vzprime=A*v-Fbar*Ghat_jth_row'; 
t(j,j)=v'*vzprime; vzprime=vzprime-t(j,j)*v;
t(j+1,j)=norm(vzprime);vzprime=vzprime/t(j+1,j); % v_2
t(j,j+1)=t(j+1,j);               
% switching location of the vectors v
v_old=v;v=vzprime;
vzprime=-v_old-W*Ghat_jth_row';
vzprime=(vzprime-t(j,j)*z)/t(j+1,j);  % z_2
% switching location of the vectors z
z_old=z;z=vzprime;


% computing angles
tau=t(1,1)'; 
c_angle_old=c_angle;          % c_0
c_angle=tau/sqrt(t(j+1,j)*t(j+1,j) + abs(tau)*abs(tau));   % c_j
s_angle=t(j+1,j)/sqrt(t(j+1,j)*t(j+1,j) + abs(tau)*abs(tau));   % s_j
gamma_old=gamma;              % gamma_0
gamma=-s_angle*gamma;         % gamma_1

% updating r, x
% r=s_angle*s_angle*r + gamma*c_angle'*v;  % r_1
x=s_angle*s_angle*x + gamma*c_angle'*z;  % x_1

j=2;

if abs(gamma)<tol,  flag=0; end; % stopping criterion for j=1

while (j<=max_it) & (flag==1)
   % j >= 2
   Fhat_jth_row=v'*[F1,G1];      % jth row of Fhat
   Ghat_jth_row=[Fhat_jth_row(s_half+1:2*s_half),-Fhat_jth_row(1:s_half)]; % jth row of Ghat 
   pstar=-s_angle*pstar + c_angle*Fhat_jth_row;    % p_j^*
   Fbar=Fbar+v*Fhat_jth_row;     % Fbar_j=f_{j:s,j} 
   W=W+z*Fhat_jth_row;           % W_j = w_{1:s,j} 
   
   
   % new Arnoldi vector
   vzprime=A*v-Fbar*Ghat_jth_row';
   t(j-1,j)=v_old'*vzprime; vzprime=vzprime-t(j-1,j)*v_old;
   t(j,j)=v'*vzprime; vzprime=vzprime-t(j,j)*v;
   t(j+1,j)=norm(vzprime);vzprime=vzprime/t(j+1,j); % v_{j+1}
   % switching location of the vectors v 
   v_old=v;v=vzprime;    
   vzprime=-v_old-W*Ghat_jth_row';
   vzprime=(vzprime-t(j,j)*z-t(j-1,j)*z_old)/t(j+1,j);  % z_{j+1}
   % switching location of the vectors z
   z_old=z;z=vzprime;
         
   % computing angles
   tau=c_angle*t(j,j) - s_angle*c_angle_old*t(j-1,j) + pstar*Ghat_jth_row'; 
   c_angle_old=c_angle;          % c_{j-1}
   c_angle=tau/sqrt(t(j+1,j)*t(j+1,j) + abs(tau)*abs(tau));   % c_j
   s_angle=t(j+1,j)/sqrt(t(j+1,j)*t(j+1,j) + abs(tau)*abs(tau));   % s_j
   gamma_old=gamma;              % gamma_{j-1}
   gamma=-s_angle*gamma;         % gamma_j

   % updating r, x
%   r=s_angle*s_angle*r + gamma*c_angle'*v;  % r_j
   x=s_angle*s_angle*x + gamma*c_angle'*z;  % x_j
   
   j=j+1;
   if abs(gamma)<tol,  flag=0; end; % stopping criterion   
end;

clear W, Fbar;
iter=j;
err=b-A*x;
