function [a,b,P,nVec] = LanTri(M,B,r)
% [a,b,P,nVec] = LanTri(M,B,r)
%
% Lanczos tridiagonalization of a complex-symmetric and block tridiagonal
% matrix. Modified partial orthogonalization is applied if required
%
% Input
%   M     M(:,:,i) are the main diagonal blocks of the block tridiagonal
%   B     B(:,:,i) are the subdiagonal blocks of the block tridiagonal
%   r     starting vector
% Outputs
%   a     main diagonal of the tridiagonal
%   b     subdiagonal of the tridiagonal
%   P     unitary
%   nVec  number of vectors selected for reorthogonalization 
% so that
%       J = P*(diag(a) + diag(b,1) + diag(b,-1))*P.'
% where J is block tridiagonal whose main diagonal blocks are M(:,:,i)
% and the subdiagonal blocks are B(:,:,i).
%
% Dependency
%    ./sbmvmul.m    symmetric block tridiagonal matrix and vector
%                   multiplication
%
n = length(r);                  % get the dimension of the starting vector
%
% constants
IM = sqrt(-1);
SQRTEPS = sqrt(eps);
MIDEPS = sqrt(SQRTEPS)^3;	% between sqrt(eps) and eps
%
% initialize two column vectors for diagonals
a = zeros(n,1); 
b = zeros(n-1,1);
wOld = zeros(n,1);		% orthogonality estimates
wCur = zeros(n,1);
wOld(1) = 1.0;
%
up = ones(n,1);	                % upper and lower bounds for
low = ones(n,1);                % orthogonalization intervals
interNum = 0;			% orthogonalization interval number
doOrtho = 0;		        % if do orthogonalization
second = 0;			% if this is the second partial orthog
nVec = 0;	                % number of vectors for reorthogonalization
%
P(:,1) = r/norm(r);             % set the first column of P
%
for j=1:n
    tmp = sbmvmul(M,B,conj(P(:,j)));    % J*conj(p(j)). Band multiplication
    a(j) = P(:,j)'*tmp;                 % a(j) = p(j)'*J*conj(p(j))
    % calculate r = J*conj(p(j)) - a(j)*p(j) - b(j-1)*p(j-1)
    if j == 1
        r = tmp - a(j)*P(:,j);
    else
        r = tmp - a(j)*P(:,j) - b(j-1)*P(:,j-1);
    end    
%
    if (j < n)
        b(j) = norm(r);     
%
        if (j > 2) 	        % compute orthogonality estimates
            wOld(1) = (b(1)*conj(wCur(2)) + a(1)*conj(wCur(1)) ...
                       - a(j)*wCur(1) - b(j-1)*wOld(1))/b(j) ...
                       + eps*(b(1)+b(j))*0.3*(randn + IM*randn);
            wOld(2:j-1) = (b(2:j-1).*conj(wCur(3:j)) ...
                           + a(2:j-1).*conj(wCur(2:j-1)) ...
                           - a(j)*wCur(2:j-1) ...
                           + b(1:j-2).*conj(wCur(1:j-2)) ...
                           - b(j-1)*wOld(2:j-1))/b(j) ...
                           + eps*0.3*(b(2:j-1)+b(j)*ones(j-2,1)) ...
                           .*(randn(j-2,1) + IM*randn(j-2,1));
%
            % swap wOld and wCur
            tmp = wOld(1:j-1);
            wOld(1:j-1) = wCur(1:j-1);
            wCur(1:j-1) = tmp;
            wOld(j) = 1.0;
        end % if j>2
        wCur(j) = eps*n*(b(1)/b(j))*0.6*(randn + IM*randn);
        wCur(j+1) = 1.0;
%
        if (second == 0)                % not the second time, determine intervals
            doOrtho = 0;	        % initialization
            interNum = 0;
            k = 1;
            while k <= j
                if (abs(wCur(k)) >= SQRTEPS)	% lost orthogonality
                    doOrtho = 1;
                    interNum = interNum + 1;
                    % find the upper bound
                    p = k + 1;
                    while ((p < (j + 1)) & (abs(wCur(p)) >= MIDEPS))
                        p = p + 1;	        % nearly lost orthogonality
                    end % while
                    up(interNum) = p - 1;
                    % find the lower bound
                    p = k - 1;
                    while ((p > 0) & (abs(wCur(p)) >= MIDEPS))
                        p = p - 1;	        % nearly lost orthogonality
                    end % while
                    low(interNum) = p + 1;
%
                    k = up(interNum) + 1;	% continue search
                else
                    k = k + 1;
                end % if lost orthogonality
            end % while k
        end % if not second time
%
        if ((doOrtho == 1) | (second == 1))	% now we have intervals,
            % carry out orthogonalization
            for (k = 1:interNum)		% for each interval
                for (i = low(k):up(k))
                    r = r - (P(:,i)'*r)*P(:,i);	% do orthogonalization
                    % reset ortho estimates
                    wCur(i) = eps*1.5*(randn + IM*randn);
                end % for i
%
                nVec = nVec + up(k) - low(k) + 1;
                % count the number of vectors selected
                if (second == 1)	        % this is the second time
                    second = 0;		        % reset
                    low(k) = 0; 
                    up(k) = 0; 
                else
                    second = 1;		        % do second time
                    doOrtho = 0;	        % reset
                    % adjust orthogonalization intervals for the second time
                    low(k) = max(1, low(k) - 1);
                    up(k) = min(j + 1, up(k) + 1);
                end % if
            end % for k
            b(j) = norm(r);                     % recalculate b(j)
        end % if
%
        if (abs(b(j)) < eps)                    % b(j)=0, quit
	    a = a(1:j); 
            b = b(1:j-1);
	    return
        else
            P(:,j+1) = r/b(j);
        end % if 
    end 
end