function [V,W,T,r,p] = NonSymmetricLanczos(A,v,n,variant)
% NONSYMMETRICLANCZOS Non-symmetric Lanczos process
%   [V,W,T,r,p]=NonSymmetricLanczos(A,v,n,variant) generates
%   biorthogonal vectors in the matrices V and W and a tridiagonal
%   matrix T with the starting vector v using the non-symmetric
%   Lanczos process, where the variant parameter can be either
%   'GolubVanLoan','Saad' or 'Gutknecht'. The auxiliary vectors r and
%   p are also returned.

V=v; W=v/norm(v)^2;
beta=0; gamma=0;
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
  switch variant
    case 'GolubVanLoan'
      beta=norm(r); gamma=r'*p/beta; 
    case 'Saad'
      beta=sqrt(abs(r'*p));gamma=beta*sign(r'*p); 
    case 'Gutknecht'
      beta=-alpha-gamma; gamma=r'*p/beta; 
  end
  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;
