www.pudn.com > NumericalComputingwithMatlabCode.zip > eigsvdgui.m, change:2003-10-28,size:7361b


function D = eigsvdgui(A,job) 
%EIGSVDGUI Demonstrate computation of matrix eigenvalues and singular values. 
%   EIGSVDGUI shows three variants of the QR algorithm. 
% 
%   EIGSVDGUI(A) for square, nonsymmetric A, or EIGSVDGUI(A,'eig'), reduces 
%   A to Hessenberg form, then applies a double-shift, eigenvalue-preserving 
%   QR algorithm.  The result is the real Schur block upper triangular form, 
%   with one-by-one diagonal blocks for real eigenvalues and two-by-two 
%   diagonal blocks for pairs of complex eigenvalues. 
% 
%   EIGSVDGUI(A) for square, symmetric A, or EIGSVDGUI(A,'symm'), reduces 
%   the symmetric part, (A+A')/2, to tridiagonal form, then applies a 
%   single-shift, eigenvalue-preserving QR algorithm.  The result is 
%   a diagonal matrix containing the eigenvalues, which are all real. 
% 
%   EIGSVDGUI(A) for rectangular A, or EIGSVDGUI(A,'svd'), reduces A to 
%   bidiagonal form, then applies a single-shift QR algorithm that preserves 
%   the singular values.  The result is a diagonal matrix containing the 
%   singular values. 
% 
%   If A is symmetric and positive definite, the three variants compute 
%   the same final diagonal matrix by three different algorithms. 
% 
%   D = EIGSVDGUI(...) returns the diagonal or Schur result. 
 
if nargin < 1 
   A = randn(24,24); 
   job = 'symm'; 
elseif nargin < 2 
   if isequal(A,A') 
      job = 'symm'; 
   elseif isequal(size(A),size(A')) 
      job = 'eig'; 
   else 
      job = 'svd'; 
   end 
elseif isequal(A,'gcf') 
   A = get(gcf,'userdata'); 
end 
 
shg 
J = jet(256); 
J(1,:) = get(gcf,'color'); 
set(gcf,'doublebuffer','on','colormap',J,'userdata',A, ... 
   'name',['eigsvdgui(A,''' job ''')'],'menu','none','numbertitle','off'); 
 
if isequal(job,'symm'); 
   A = eiggui((A+A')/2); 
elseif isequal(job,'eig') 
   A = eiggui(A); 
else 
   A = svdgui(A); 
end 
 
eig = uicontrol('units','norm','pos',[.02,.02,.10,.04], ... 
   'string','eig','callback','eigsvdgui(''gcf'',''eig'')'); 
symm = uicontrol('units','norm','pos',[.14,.02,.10,.04], ... 
   'string','symm', 'callback','eigsvdgui(''gcf'',''symm'')'); 
svd = uicontrol('units','norm','pos',[.26,.02,.10,.04], ... 
   'string','svd', 'callback','eigsvdgui(''gcf'',''svd'')'); 
stop = uicontrol('units','norm','pos',[.38,.02,.10,.04], ... 
   'string','close','callback','close'); 
if ~isequal(size(A),size(A')) 
   set([eig,symm],'foreground',[.66 .66 .66],'callback',[]) 
end 
if nargout > 0 
   D = A; 
end 
 
 
% ------------------------------------------- 
 
function A = eiggui(A) 
 
scale = 256/sqrt(max(abs(diag(A'*A)))); 
imageh = image(ceil(scale*abs(A))+1); 
daspect([1 1 1]) 
issymm = isequal(A,A'); 
iscmplx = ~isreal(A); 
 
% Househoulder reduction to tridiagonal or Hessenberg form. 
 
[n,n] = size(A); 
for k = 1:n-2 
 
   % Introduce zeros below the subdiagonal in the k-th column. 
 
   u = A(:,k); 
   u(1:k) = 0; 
   sigma = norm(u); 
   if sigma ~= 0 
      if u(k+1) ~= 0, sigma = sign(u(k+1))*sigma; end 
      u(k+1) = u(k+1) + sigma; 
      rho = 1/(sigma'*u(k+1)); 
      v = rho'*A*u; 
      w = (rho*u'*A)'; 
      gamma = rho/2*u'*v; 
      v = v - gamma*u; 
      gamma = rho/2*u'*w; 
      w = w - gamma*u; 
      A = A - v*u' - u*w'; 
      A(k+2:n,k) = 0; 
      if issymm, A(k,k+2:n) = 0; end 
   end 
   set(imageh,'cdata',ceil(scale*abs(A))+1) 
   pause(.1) 
end 
 
% Tridiagonal or Hessenberg QR algorithm. 
 
it = 0; 
titleh = title('0'); 
k = n; 
while k > 1 
 
   % 1-by-1 convergence test. 
 
   if abs(A(k,k-1)) <= 2*eps*(abs(A(k-1,k-1)) + abs(A(k,k))) 
      A(k,k-1) = 0; 
      if issymm 
         A(k-1,k) = 0; 
         A(k,k) = real(A(k,k)); 
      end 
      k = k-1; 
   else 
 
      % Wilkinson shift, eigenvalues of lower 2-by-2, A(k-1:k,k-1:k). 
    
      r = (A(k,k)-A(k-1,k-1))/(2*A(k,k-1)); 
      s = r^2 + A(k-1,k)/A(k,k-1); 
    
      % Use single shift for real eigenvalues of real matrices 
      % and for all eigenvalues of complex matrices. 
    
      if iscmplx | s >= 0  
    
         % Single real shift, eigenvalue of 2-by-2 closest to A(k,k). 
    
         s = sqrt(s); 
         if r < 0, s = -s; end 
         if r+s ~= 0, s = A(k,k) + A(k-1,k)/(r+s); end 
    
         % Single QR step. 
    
         I = eye(k,k); 
         [Q,R] = qr(A(1:k,1:k) - s*I); 
         A(1:k,1:k) = R*Q + s*I; 
         it = it+1; 
    
      else 
    
         % Complex eigenvalues of real matrices. 
         % 2-by-2 convergence test. 
    
         if k == 2 
            k = 0; 
         elseif abs(A(k-1,k-2)) <= 2*eps*(abs(A(k-2,k-2)) + abs(A(k-1,k-1))) 
            A(k-1,k-2) = 0; 
            if issymm, A(k-2,k-1) = 0; end 
            k = k-2; 
         else 
    
            % Sum and product of eigenvalues of lower 2-by-2. 
       
            t = A(k-1,k-1) + A(k,k); 
            d = A(k-1,k-1)*A(k,k) - A(k,k-1)*A(k-1,k); 
       
            % Double QR step. 
       
            I = eye(k,k); 
            [Q,R] = qr(A(1:k,1:k)^2 - t*A(1:k,1:k) + d*I); 
            A(1:k,1:k) = triu(Q'*A(1:k,1:k)*Q,-1); 
            it = it+2; 
         end 
      end 
   end 
   if issymm, A(1:k,1:k) = tril(A(1:k,1:k),1); end 
   set(imageh,'cdata',ceil(scale*abs(A))+1) 
   set(titleh,'string',num2str(it)) 
   pause(.1) 
end 
if issymm, A(1,1) = real(A(1,1)); end 
 
 
% ------------------------------------------- 
 
function A = svdgui(A) 
%SVDGUI Demonstrate the computation of the SVD. 
%   SVDGUI(A) shows the steps in the computation of the 
%   singular value decomposition of any real or complex matrix. 
 
scale = 256/sqrt(max(abs(diag(A'*A)))); 
imageh = image(ceil(scale*abs(A))+1); 
daspect([1 1 1]) 
 
% Househoulder reduction to bidiagonal form. 
 
[m,n] = size(A); 
for k = 1:min(m,n) 
 
   % Introduce zeros below the diagonal in the k-th column. 
 
   u = A(:,k); 
   u(1:k-1) = 0; 
   sigma = norm(u); 
   if sigma ~= 0 
      if u(k) ~= 0, sigma = sign(u(k))*sigma; end 
      u(k) = u(k) + sigma; 
      rho = 1/(sigma'*u(k)); 
      v = rho*(u'*A); 
      A = A - u*v; 
      A(k+1:m,k) = 0; 
   end 
   set(imageh,'cdata',ceil(scale*abs(A))+1); 
   pause(.1) 
 
   % Introduce zeros to the right of the superdiagonal in the k-th row. 
 
   u = A(k,:); 
   u(1:k) = 0; 
   sigma = norm(u); 
   if sigma ~= 0 
      if u(k+1) ~= 0, sigma = sign(u(k+1))*sigma; end 
      u(k+1) = u(k+1) + sigma; 
      rho = 1/(sigma'*u(k+1)); 
      v = rho*(A*u'); 
      A = A - v*u; 
      A(k,k+2:n) = 0; 
   end 
   set(imageh,'cdata',ceil(scale*abs(A))+1); 
   pause(.1) 
end 
 
% Bidiagonal SVD QR iteration. 
 
it = 0; 
titleh = title('0'); 
k = min(m,n); 
while k > 1 
 
   % Convergence test. 
 
   if abs(A(k-1,k)) <= 2*eps*(abs(A(k-1,k-1)) + abs(A(k,k))) 
      A(k-1,k) = 0; 
      k = k-1; 
   else 
 
      % One step of single shift QR iteration. 
      % Wilkinson shift, eigenvalue of lower 2-by-2 of A'*A. 
    
      T = A(1:k,1:k)'*A(1:k,1:k); 
      r = (T(k,k)-T(k-1,k-1))/(2*T(k,k-1)); 
      s = sqrt(r^2 + T(k-1,k)/T(k,k-1)); 
      if r < 0, s = -s; end 
      if r+s ~= 0, s = T(k,k) + T(k-1,k)/(r+s); end 
 
      I = eye(k,k); 
      [Q,R] = qr(T-s*I); 
      A(1:k,1:k) = A(1:k,1:k)*Q; 
      [Q,R] = qr(A(1:k,1:k)); 
      A(1:k,1:k) = tril(R,1); 
      it = it+1; 
   end 
 
   set(imageh,'cdata',ceil(scale*abs(A))+1); 
   set(titleh,'string',num2str(it)) 
   pause(.1) 
end 
A = abs(A);