function [V,W,T,f,g] = NSLanczos(A,k,b,c);
%
%   Usage:  [V,W,T,f,g] = NSLanczos(A,k,b,c);
%
%   Input:  A -- an n by n matrix  
%           k -- a positive integer (k << n assumed)
%           b -- an n  vector (b'*c .ne. 0 assumed)
%           c -- an n  vector (b'*c .ne. 0 assumed)
%
%   Output: V -- an n by k bi-orthogonal matrix
%           W -- an n by k bi-orthogonal matrix  with W'*V = I_k
%           T -- a k by k  non-symmetric tridiagonal matrix
%           f -- an n vector  (W'*f = 0)
%           g -- an n vector  (V'*g = 0)
%
% 
%           with   AV = VT + fe_k', and   A'W = WT' + ge_k'
%
%   In real life you would not store V, and you would store T as three 
%   vectors  a = diag(T) , b = diag(T,-1), c = diag(T,1);
%
%   D.C. Sorensen
%   05 Nov 03
%
    n = length(b);
    T = zeros(k);
    V = zeros(n,k); W = V;

    v1 = b/sqrt(abs(b'*c));
    scl = v1'*c;
    w1 = c/scl;
 
    f = A*v1;
    alpha = w1'*f;
    f = f - v1*alpha;  g = A'*w1 - w1*alpha;

    V(:,1) = v1;  W(:,1) = w1;  T(1,1) = alpha;

    for j = 2:k,

        gamma = g'*f;  

        if (abs(gamma) == 0), disp('breakdown'); return; end

        beta = sqrt(abs(gamma)); gamma = gamma/beta;

        v0 = v1; v1 = f/beta;    w0 = w1; w1 = g/gamma;

        f = A*v1 - v0*gamma;   

        alpha = w1'*f;

        f = f - v1*alpha;  
        
        g = (A'*w1 - w0*beta) - w1*alpha;

        T(j,j-1) = beta; T(j-1,j) = gamma; T(j,j) = alpha;
        V(:,j)   = v1; W(:,j)   = w1;

    end