function [d,iters,flopi,Z] = imtql(d,e,Z)
%   [d,iters,flopi] = imtql(d,e)  % No Z
%   [d,iters,flopi,Z] = imtql(d,e)  % initial Z = I
%   [d,iters,flopi,Z] = imtql(d,e,Z)

%   Eispack subroutine imtql2(nm,n,d,e,z,ierr)
%   This subroutine is a translation of the Algol procedure imtql2,
%   Num. Math. 12, 377-383(1968) by Martin and Wilkinson,
%   as modified in Num. Math. 15, 450(1970) by Dubrulle.
%   Handbook for Auto. Comp., vol.ii-Linear Algebra, 241-248(1971).

    n = length(d);
    if length(e) == n
        e(1:n-1) = e(2:n);
    end
    e(n) = 0;
    
    wantZ = nargout > 3;
    if wantZ && nargin < 3
        Z = eye(n,n);
    end
   
    iters = zeros(size(d));
    flopi = zeros(size(d));
    for k = 1:n,                                                flops("0")
        while 1  % iteration
            m = k;
            while m < n
                % look for small sub-diagonal element
                tst1 = abs(d(m)) + abs(d(m+1)); 
                tst2 = tst1 + abs(e(m));                        flops(2)
                if tst2 == tst1
                    break
                end
                m = m + 1;
                p = d(k);
            end
            if m == k
                % exit iteration loop
                break
            end
            if iters(k) == 30
                error(['**** 30 iterations for k = ' int2str(k)])
            end
            iters(k) = iters(k) + 1;
            % form shift
            g = (d(k+1) - p) / (2*e(k));
            r = hypot(g,1.0);
            g = d(m) - p + e(k) / (g + sign(g+realmin)*r);      flops(8)
            % chase bulge
            s = 1.0;
            c = 1.0;
            p = 0.0;
            for i = m-1:-1:k
                f = s * e(i);
                b = c * e(i);
                r = hypot(f,g);                                 flops(3)
                e(i+1) = r;
                if r == 0.0
                    % underflow
                    d(i+1) = d(i+1) - p;
                    e(m) = 0.00;
                    break
                end
                s = f / r;
                c = g / r;
                g = d(i+1) - p;
                r = (d(i) - g) * s + 2.0 * c * b;
                p = s * r;
                d(i+1) = g + p;  
                g = c * r - b;                                 flops(10)
                if wantZ
                    % eigenvectors
                    Z(:,i:i+1) = Z(:,i:i+1) * [c s; -s c];     flops(4*n)
                end
            end
            d(k) = d(k) - p;                                   flops(1)
            e(k) = g;
            e(m) = 0;
        end % iteration
                                                    flopi(k) = flops;
    end % for k
    
    [d,p] = sort(d);
    iters = iters(p);
    flopi = flopi(p);
    if wantZ
        Z = Z(:,p);
    end

    function fl = flops(a)
        persistent count
        if nargin == 0
            fl = count;
        elseif isstring(a)
            count = 0;
        else
            count = count + a;
        end
    end
end