% QMR algorithm 
% bi-orthogonal Lanczos procedure
%
% A(N,N), f(N,1), x0(N,1) given then 
% 
% v(i+1)= A v(i) - a(i) v(i) - b(i) v(i-1)
% w(i+1)= A'w(i) - a'(i)w(i) - c(i) w(i-1)
% 
% with a, b, c such that w(i)'*v(i+1)=w(i-1)'*v(i+1)=0
%
% then normalization: 
%  d=norm(v(i+1)); v(i+1)=v(i+1)/d;   c(i+1)=d ;
%  e=(w1'*v1);     w(i+1)=w(i+1)/e';  b(i+1)=e ;
% 
% F. Waleffe 11/30/2000

% Initialization

x=x0;  r0=f-A*x; rnrm=norm(r0); resvec=[rnrm]; N=length(x); maxit=N;

v1=r0/rnrm; w1=conj(v1); g=[rnrm;0];

b=0; c=0; v0=zeros(N,1); w0=v0; p0=v0; p1=p0;  % init v, w, p 3-term recurrence
G0=eye(2); G1=G0;                              % init to no rotation

for i=1:maxit,

   v2=A*v1; w2=A'*w1; a=w1'*v2; 
   
   v2=v2-a*v1-b*v0; w2=w2-a'*w1-c*w0;         % bi-orthogonalization
   c=norm(v2); v2=v2/c;                       % v normalization
   
   R=[0;b;a;c]; 
   R(1:2)=G0*R(1:2);     R(2:3)=G1*R(2:3);    % apply old rotations
   G2=givens(R(3),R(4));                      % new rotation
   R(3:4)=G2*R(3:4);     g=G2*g;              % apply new rotation
   
   p2=(v1-R(2)*p1-R(1)*p0)/R(3);              % search directions
   x = x + g(1)*p2;         
   res=norm(f-A*x);  resvec=[resvec;res];
      
   b=(w2'*v2); w2=w2/b';                      %  w normalization
   
   v0=v1; w0=w1; v1=v2; w1=w2;                % update
   G0=G1; G1=G2; p0=p1; p1=p2; g=[g(2);0];

end

