function [a,b,V,U] = gLanPO(A,rv)

% [a,b,V,U] =gLanPO(A,rv)

%

% General Lanczos bidiagonalization with partial orthogonalization.

% Given a (complex) matrix A, this computes

% a, b, and unitary V and U such that

%       A*V = U*(diag(a) + diag(b,1))

%

% Input

%   A	m-by-n (m >= n) complex matrix

%   rv  n-by-1, start vector for V

% Outputs

%   a,b diagonal and subdiagonal

%   V   n-by-n unitary matrix

%   U   m-by-n with orthnormal columns



% S. Qiao	McMaster Univ. May 2007 

%

% Reference

% Gene H. Golub and Charles F. Van Loan.

% Matrix Computations, 3rd Edition.

% The Johns Hopkins University Press,

% Baltimore, MD. 1996. p. 495.



[m,n] = size(A);



% constants

IM = sqrt(-1);

SQRTEPS = sqrt(eps);



a = zeros(n,1);

b = zeros(n-1,1);

V = zeros(n,n);

U = zeros(m,n);



p = rv/norm(rv);	% initialize a unit 2-norm vector

bOld = 1.0; uOld = zeros(m,1);



nuOld = zeros(n,1);     % initial orthogonality estimates

nuOld(1) = 1.0;

nuCur = zeros(n,1);

muOld = zeros(n,1);

muCur = zeros(n,1);

muCur(1) = eps*m*0.6*(randn + IM*randn);

muCur(2) = 1.0;



secondU = 0; secondV = 0; % if do ortho in the next iteration

nVecV = 0; nVecU = 0;   % no of vectors selected for ortho

%

for j=1:n

    vj = p/bOld;

    V(:,j) = vj;

% find a column of U

    r = A*vj - bOld*uOld;

    a(j) = norm(r);

% estimate the orthogonality of r (U(:,j)) against U(:,1:j-1)

    if (j > 2)

        muOld(1:j-2) = (conj(b(1:j-2)).*nuCur(2:j-1) ...

                      + conj(a(1:j-2)).*nuCur(1:j-2) ...

                      - b(j-1)*muCur(1:j-2))/a(j) ...

       + eps*(b(1:j-2) + a(j)*ones(j-2,1))*0.3*(randn + IM*randn);

        % swap muOld(1:j-2) and muCur(1:j-2)

        tmp = muOld(1:j-2);

        muOld(1:j-2) = muCur(1:j-2);

        muCur(1:j-2) = tmp;

        muOld(j-1) = 1.0;

        muCur(j-1) = eps*m*((a(2)*b(1)')/(a(j)*b(j-1)')) ...

                    *0.6*(randn + IM*randn);

        muCur(j) = 1.0;

    end % if j>2



    if ((max(abs(muCur(1:j-1))) > SQRTEPS) | (secondU==1))

                               % orthogonalization

        if (max(abs(muCur(1:j-1))) > SQRTEPS)

                               % if orthogonality is lost

            secondU = 1;       % do it in the next iteration

        else                   % if this the second time

            secondU = 0;       % reset

        end



        for k = 1:j-1          % MGS orthogonalization

            r = r - (U(:,k)'*r)*U(:,k);

            muCur(k) = eps*1.5*(randn + IM*randn);

                               % adjust estimates

        end % for k



        nVecU = nVecU + j - 1;

        a(j) = norm(r);    % recalculate a(j)

    end % if ((max ...



    uj = r/a(j);

    uOld = uj;

    U(:,j) = uj;



    % find a column of V

    if (j < n)

        p = A'*uj - a(j)*vj;

        b(j) = norm(p);

% estimate the orthogonality of p (V(:,j+1)) against V(:,1:j)

        if (j > 2)

            nuOld(1) = (a(1)'*muCur(1) - a(j)*nuCur(1))/b(j) ...

                  + eps*(a(1) + b(j))*0.3*(randn + IM*randn);

            nuOld(2:j-1) = (conj(a(2:j-1)).*muCur(2:j-1) ...

                          + conj(b(1:j-2)).*muCur(1:j-2) ...

                          - a(j)*nuCur(2:j-1))/b(j) ...

       + eps*(a(2:j-1) + b(j)*ones(j-2,1))*0.3*(randn + IM*randn);

            % swap nuOld(1:j-1) and nuCur(1:j-1)

            tmp = nuOld(1:j-1);

            nuOld(1:j-1) = nuCur(1:j-1);

            nuCur(1:j-1) = tmp;

            nuOld(j) = 1.0;

        end % if j>2

        nuCur(j) = eps*n*((a(1)*b(1)')/(a(j)*b(j)')) ...

                        *0.6*(randn + IM*randn);

        nuCur(j+1) = 1.0;



        if ((max(abs(nuCur(1:j))) > SQRTEPS) | (secondV==1))

                               % orthogonalization

            if (max(abs(nuCur(1:j))) > SQRTEPS)

                               % if orthogonality is lost

                secondV = 1;   % do it in the next iteration

            else               % if this the second time

                secondV = 0;   % reset

            end



            for k = 1:j        % MGS orthogonalization

                p = p - (V(:,k)'*p)*V(:,k);

                nuCur(k) = eps*1.5*(randn + IM*randn);

                               % adjust estimates

            end % for k



            nVecV = nVecV + j;

            b(j) = norm(p);    % recalculate b(j)

        end % if ((max ...

        bOld = b(j);

    end % if j<n



end % for j

%fprintf('\n %d U vectors selected for orthogonalization.', nVecU)

%fprintf('\n %d V vectors selected for orthogonalization.', nVecV)