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

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

%

% General Lanczos bidiagonalization with modified partial orthogonalization 

% and restart.

% 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.



% is A complex?

if (imag(A) ~= 0)

    isComplex = 1;

else

    isComplex = 0;

end



[m,n] = size(A);



% constants

IM = sqrt(-1);

SQRTEPS = sqrt(eps);

MIDEPS = sqrt(SQRTEPS)^3; % between sqrt(eps) and eps

RSTOL = SQRTEPS*(norm(A,'fro')/(m*n)); % tolerance for restart

MAXRS = 50; % maximal number of restarts



a = zeros(n,1);

b = zeros(n-1,1);

V = zeros(n,n);

U = zeros(m,n);





if (norm(rv) == 0)

    rv = ones(n,1) - 2*rand(n,1);

    if (isComplex == 1)

        rv = rv + IM*(ones(n,1) - 2*rand(n,1));

    end

end

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



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

normp = 1.0;



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;    % number of vectors selected for orthog

upV(1) = 0; lowV(1) = 0; upU(1) = 0; lowU(1) = 0;

%

for j=1:n



    vj = p / normp;

    V(:,j) = vj;

    

    r = A*vj - bOld*uOld;    % find a column of U



    a(j) = norm(r);



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

    if (j > 2)

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

                     + 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 (jc > 2)



    % not the second time, determine intervals

    if (secondU == 0)	

        interNumU = 0;

        k = 1;

        while k <= j

            if (abs(muCur(k)) > SQRTEPS)    % lost orthogonality

                interNumU = interNumU + 1;

                % find the upper bound

                curpos = k + 1;

                % nearly lost orthogonality

                while ((curpos <= (j - 1)) & (abs(muCur(curpos)) >= MIDEPS))

                    curpos = curpos + 1;    

                end % while

                upU(interNumU) = curpos - 1;

                % find the lower bound

                curpos = k - 1;

                % nearly lost orthogonality

                while ((curpos > 0) & (abs(muCur(curpos)) >= MIDEPS))

                    curpos = curpos - 1;    

                end % while

                lowU(interNumU) = curpos + 1;

                %

                k = upU(interNumU) + 1;    % continue search

            else

                k = k + 1;

            end % if lost orthogonality

        end % while k

    end % if not second time



    % now we have intervals, carry out orthogonalization

    if (interNumU > 0 | (secondU == 1))

        for (k = 1:interNumU)    % for each interval

	    for (i = lowU(k):upU(k))    % do orthogonalization

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

                muCur(i) = eps*1.5*(randn + IM*randn);    % reset ortho estimates

            end % for i



            nVecU = nVecU + upU(k) - lowU(k) + 1;    % count the vectors selected

        

            if (secondU == 1)    % this is the second time

                lowU(k) = 0; upU(k) = 0; 

            else    % this is the first time

                % adjust orthog intervals for the second time

                lowU(k) = max(1, lowU(k) - 1);

                upU(k) = min(j, upU(k) + 1);

            end % if second == 1

        end % for k

    

        secondU = ~secondU;    % repeat if not already second time

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

    end % if interNum > 0 ...



    uj = r/a(j);

    uOld = uj;

    U(:,j) = uj;

    



    % find a column of V

    if (j < n)

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



        normp = norm(p);

        b(j) = normp;

        

        % determine orthogonalization intervals

	nRS = 0;

	nEntry = 0;    % number of entries into the following loop

	while ((nEntry < 1) | ((normp < RSTOL) & (nRS < MAXRS)))

	    nEntry = nEntry + 1;

	    if (normp < (RSTOL * n * m))  % small subdiagonal detected

%                fprintf('\nSmall subdiagonal encountered. Restart V ...  (%d   %d)',normp, (RSTOL * n * m));

	        % orthogonalize p against all previous V(:,1:j)

	        interNumV = 1;    % one ortho interval [1:j]

	        lowV(interNumV) = 1;

	        upV(interNumV) = j;

	        p = ones(n,1) - 2*rand(n,1);    % restart

                if (isComplex == 1)

                    p = p + IM*(ones(n,1) - 2*rand(n,1));

                end

	        nRS = nRS + 1;

	        % update nuOld and nuCur

	        if (j >=2)

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

                    nuOld(j) = 1.0;

                end

                nuCur(j+1) = 1.0;

            else % subdiagonal is not small, no restart               



                % 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) = (a(2:j-1).*muCur(2:j-1) ...

                               + 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 (secondV == 0)    % not the second time, determine intervals

                    interNumV = 0;

                    k = 1;

                    while k <= (j-1)

                        if (abs(nuCur(k)) > SQRTEPS)    % lost orthogonality

                            interNumV = interNumV + 1;

                            % find the upper bound

                            curpos = k + 1;

                            % nearly lost orthogonality

                            while ((curpos < (j + 1)) & ...

                                   (abs(nuCur(curpos)) >= MIDEPS))

                                curpos = curpos + 1;    

                            end % while

                            upV(interNumV) = curpos - 1;

                            % find the lower bound

                            curpos = k - 1;

                            % nearly lost orthogonality

                            while ((curpos > 0) & ...

                                   (abs(nuCur(curpos)) >= MIDEPS))

                                curpos = curpos - 1;    

                            end % while

                            lowV(interNumV) = curpos + 1;

                            %

                            k = upV(interNumV) + 1;    % continue search

                        else

                            k = k + 1;

                        end % if lost orthogonality

                    end % while k

                end % if not second time

            end



            % now we have intervals, carry out orthogonalization

            if (interNumV > 0 | (secondV == 1))	

                % for each interval, do orthog, reset ortho estimates

                for (k = 1:interNumV) 

                    for (i = lowV(k):upV(k))

                        p = p - (V(:,i)'*p)*V(:,i); % do orthogonalization

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

                    end % for i



                    nVecV = nVecV + upV(k) - lowV(k) + 1;

                    

                    % count the number of vectors selected

                    if (secondV == 1)	% this is the second time

                        lowV(k) = 0; upV(k) = 0; 

                    else

                        % adjust orthog intervals for the second time

                        lowV(k) = max(1, lowV(k) - 1);

                        upV(k) = min(j + 1, upV(k) + 1);

                    end % if

                end % for k

                secondV = ~secondV;    % repeat if not already second time	                    

                normp = norm(p);

            end % if    

        end        

        

        bOld = b(j);

    

        if (normp < RSTOL)

            fprintf('\nRestart %d times', nRSV);

	    a = a(1:j); 

	    b = b(1:j-1);

	    nsteps = j;

	    return

	end

	

	% if no restart, we may have performed orthogonalization,

	% recalculate b(j); if restart, we use the old small b(j)

	% instead of normr.

	if (nRS == 0)       % no restart

	    b(j) = normp;   % recalculate b(j)

        end

       

    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)