function [X,res,V,W,T]=Saad(A,b,x0,n);
% SAAD Non-symmetric Lanczos for linear systems
%   [X,res,V,W,T]=Saad(A,b,x0,n) computes an approximate solution of
%   the linear system Ax=b performing n non-symmetric Lanczos
%   iterations starting with the initial vector x0. The matrix X
%   contains the approximate solutions and the vector res the norms of
%   the corresponding residuals. The method also returns the
%   biorthogonal vectors in V and W and the tridiagonal matrix T

X=x0; r0=b-A*x0; res=norm(r0);
V=r0; W=r0/norm(r0)^2;
beta=0; gamma=0; r=r0;
alpha=W(:,1)'*A*V(:,1);
T(1,1)=alpha;
for j=1:n,
  z=T\eye(j,1); 
  x=x0+V*z;
  X=[X x];
  if j==1,
    r=A*V(:,j)-alpha*V(:,j);
    p=A'*W(:,j)-alpha*W(:,j);
  else
    r=A*V(:,j)-alpha*V(:,j)-gamma*V(:,j-1);
    p=A'*W(:,j)-alpha*W(:,j)-beta*W(:,j-1);
  end
  beta=sqrt(abs(r'*p)); gamma=beta*sign(r'*p); 
  res=[res norm(r)*abs(z(j))];  
  T(j+1,j)=beta; T(j,j+1)=gamma;
  V=[V, r/beta]; W=[W, p/gamma];
  alpha=(W(:,j+1)'*A*V(:,j+1));
  T(j+1,j+1)=alpha;
end;