函数[lambda，erers，flopz，q] = symeig（a）％symeig使用wilkinson的tred和imtql来计算％一个真正symmetrix矩阵的特征值，a。％[lambda，erters，flopz，q] = symeig（a）％lambda，特征值％erers，迭代计数％flopz，浮点操作计数％q，正交矩阵，q'* a * q = diag（lambda），q'* q = i％q'* a * q = diag（d）+诊断（e，1）+诊断（e，-1）[d，e，flopt，q] = tred（a）;[d，erats，flopi，q] = imtql（d，e，q）;flopz = flopt + sum（flopi）;erats = sum（erters）;lambda = D（:);％-----------------------------------------------------------功能[d，e，flopt，z] = try（a）％[d，e，flopt，z] = tred2（a）%% eispack子程序TRED2（NM，N，A，D，E，Z）％该子程序是藻类程序TRED2，％NUM的转换。数学。 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