function [lambda,iters,flopz,Q] = symeig(A)
% symeig Use Wilkinson's tred and imtql to compute the eigenvalues of
% a real symmetrix matrix, A.
% [lambda,iters,flopz,Q] = symeig(A)
% lambda, eigenvalues
% iters, iteration count
% flopz, floating point operations count
% Q, orthgonal matrix, Q'*A*Q = diag(lambda), Q'*Q = I
% Q'*A*Q = diag(d) + diag(e,1) + diag(e,-1)
[d,e,flopt,Q] = tred(A);
[d,iters,flopi,Q] = imtql(d,e,Q);
flopz = flopt + sum(flopi);
iters = sum(iters);
lambda = d(:);
% -----------------------------------------------------------------
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;
end
function [d,iters,flopi,Z] = imtql(d,e,Z)
% [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);
e(1:n-1) = e(2:n);
e(n) = 0;
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)
s = 1.0;
c = 1.0;
p = 0.0;
% chase bulge
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)
% eigenvectors
Z(:,i:i+1) = Z(:,i:i+1) * [c s; -s c]; flops(4*n)
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);
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