function [d,e,flopt,Z] = tred(A)
%   [d,e,flopt,Z] = tred2(A)
%
%   Eispack subroutine tred2(nm,n,a,d,e,z)
%   This subroutine is a translation of the Algol procedure tred2,
%   Num. Math. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson,
%   Handbook for Auto. Comp., vol.ii-Linear Algebra, 212-226(1971). 

   [~,n] = size(A);
   Z = A;
   d = A(n,:);
   e = zeros(size(d));                                 flops("0")
   for i = n:-1:2
       h = 0;
       scale = sum(abs(d(1:i-1)));
       if scale == 0
           e(i) = d(i-1);
           j = 1:i-1;
           Z(j,i) = 0;
       else
           k = 1:i-1;
           d(k) = d(k)/scale;                          flops(i-1)
           h = sum(d(k).^2);                           flops(i)              
           f = d(i-1);
           g = -sign(f)*sqrt(h);     
           e(i) = scale*g;                      
           h = h - f*g;                   
           d(i-1) = f - g;                             flops(5)
           e(1:i-1) = 0;
           for j = k
               f = d(j);
               Z(j,i) = f;
               g = e(j) + Z(j,j)*f;                    flops(2)
               for k = j+1:i-1
                   g = g + Z(k,j)*d(k);      
                   e(k) = e(k) + Z(k,j)*f;             flops(4)
               end
               e(j) = g;
           end
       end
       j = 1:i-1;
       e(j) = e(j)/h;                                  flops(i-1)
       f =  e(j)*d(j)';                                flops(i-1)
       hh = f/(h + h);                                 flops(1)
       e(j) = e(j) - hh*d(j);                          flops(2*i-2)
       for j = 1:i-1
           k = j:i-1;
           Z(k,j) = Z(k,j) - d(j)*e(k)' - e(j)*d(k)';  flops(4*(i-j))
           d(j) = Z(i-1,j);
           Z(i,j) = 0;
       end
       d(i) = h;
   end
   for i = 2:n
       Z(n,i-1) = Z(i-1,i-1);
       Z(i-1,i-1) = 1;
       h = d(i);
       k = (1:i-1)';
       if h ~= 0
           d(k) = Z(k,i)/h;                            flops(1)
           for j = k
               g = Z(k,i)'*Z(k,j);                     flops(i-1)
               Z(k,j) = Z(k,j) - g.*d(k)';             flops(2*i-2)
           end
       end
       Z(k,i) = 0;
   end
   d = Z(n,:);
   Z(n,:) = 0;
   Z(n,n) = 1;
   e(1) = 0;                                   flopt = flops;

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