function [X,res,V,W,T]=BiORES(A,b,x0,n);
% BIORES Non-symmetric Lanczos for linear equations
%   [X,res,V,W,T]=BiORES(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. X contains the
%   successive approximations computed by the proposal of Gutknecht,
%   and the vector res contains 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,
  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=-alpha-gamma; 
  if j==1
    x=-(V(:,j)+alpha*X(:,j))/beta;
  else
    x=-(V(:,j)+alpha*X(:,j)+gamma*X(:,j-1))/beta;
  end
  X=[X x];
  gamma=r'*p/beta;
  res=[res norm(r/beta)];
  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