function m = blocking(A,delta,showplot)
%BLOCKING  Produce blocking pattern for block Parlett recurrence.
%          M = BLOCKING(A, DELTA, SHOWPLOT) accepts an upper triangular matrix
%          A and produces a blocking pattern, specified by the vector M,
%          for the block Parlett recurrence.
%          M(i) is the index of the block into which A(i,i) should be placed.
%          DELTA is a gap parameter (default 0.1) used to determine the
%          blocking.
%          Setting SHOWPLOT nonzero produces a plot of the eigenvalues
%          that indicates the blocking:
%            - Black circles show a set of 1 eigenvalue.
%            - Blue circles show a set of >1 eigenvalues.
%              The lines connect eigenvalues in the same set.
%              Red squares show the mean of each set.

%          For A coming from a real matrix it should be posible to take
%          advantage of the symmetry about the real axis.  This code does not.

a = diag(A); n = length(a);
m = zeros(1,n); maxM = 0;

if nargin < 2 | isempty(delta), delta = 0.1; end
if nargin < 3, showplot = 0;  end

if showplot, clf, hold on, end

for i = 1:n

    if m(i) == 0
        m(i) = maxM + 1; % If a(i) hasn`t been assigned to a set
        maxM = maxM + 1; % then make a new set and assign a(i) to it.
    end

    for j = i+1:n
        if m(i) ~= m(j)    % If a(i) and a(j) are not in same set.
            if abs(a(i)-a(j)) <= delta
                if showplot
                    plot(real([a(i) a(j)]),imag([a(i) a(j)]),'o-')
                end

                if m(j) == 0
                    m(j) = m(i); % If a(j) hasn`t been assigned to a
                                 % set, assign it to the same set as a(i).
                else
                    p = max(m(i),m(j)); q = min(m(i),m(j));
                    m(m==p) = q; % If a(j) has been assigned to a set
                                 % place all the elements in the set
                                 % containing a(j) into the set
                                 % containing a(i) (or vice versa).
                    m(m>p) = m(m>p) -1;
                    maxM = maxM - 1;
                                 % Tidying up. As we have deleted set
                                 % p we reduce the index of the sets
                                 % > p by 1.
                end
            end
        end
    end
end

if showplot
    for i = 1:max(m)
        a_ind = a(m==i);
        if length(a_ind) == 1
            plot(real(a_ind),imag(a_ind),'ok' )
        else
%            plot(real(mean(a_ind)),imag(mean(a_ind)),'sr' )
        end
    end
    grid
    hold off
    box on
end