powereig.m

matlab算法集 matlab算法集

function [c,x,k] = powereig (A,tol,m,dir,dm)
%----------------------------------------------------------------
% Usage:       [c,x,k] = powereig (A,tol,m,dir,dm)
%
% Description: Apply the power method or the inverse power method
%              to find the dominant eigenvalue of a matrix and 
%              its eigenvector or the least dominant eigenvalue 
%              and its eigenvector.
%
% Inputs:      A   = n by n coefficient matrix 
%              tol = error tolerance used to terminate the
%                    search (tol >= 0).
%              m   = maximum number of iterations (m >= 1)
%              dir = direction code defined as follows:
%
%                     1 = find dominant eigenvalue
%                    -1 = find least dominant eigenvalue
%
%              dm  = optional display mode.  If present,
%                    intermediate results are displayed.
%
% Outputs:     c   = eigenvalue: det(cI - A) = 0
%              x   = n by 1 eigenvector: Ax = cx
%              k   = number of iterations performed.  If 0<k<m,
%                    then the following termination criterion
%                    was satisfied where r(x) = (A - cI)x is
%                    the error vector:
%
%                    ||r(x)|| < tol
%
% Notes:       For powereig to be successful when dir = 1, A must
%              have a single dominant eigenvalue. When dir = -1,
%              A must be nonsingular and have a single least
%              dominant eigenvalue.  If A dir = -1 and A is 
%              singular, then k is set to zero.
%----------------------------------------------------------------

% Initialize

   tol = args (tol,0,tol,2,'powereig');
   m   = args (m,1,m,3,'powereig');
   display = nargin > 4;
   n = size (A,1);
   r = zeros (n,1);
   y = randu (n,1,-1,1);
   k = 0;
   t = tol + 1;

   if dir >= 0 

% Power method 

      x = A*y;
      while (t > tol) & (k < m) 
         y = x/norm(x,inf);
         x = A*y;
         c = dot(y,x)/dot(y,y);
         r = c*y - x;
         t = norm (r,inf);
         k = k + 1; 
         if display
            fprintf ('\n%g & %.7f',k,c);
            fprintf (' & %.7f',y);
            fprintf (' & %.7f \\\\',t);
         end
      end
      
   else 

% Inverse power method 

      if abs(det(A)) < eps                    % singular A
         k = 0;
         c = 0;
         fprintf ('\nThe matrix in powereig is singular.\n');
         wait 
         return;
      else                                    % nonsingular A
         Q = zeros (n,n);                     
         p = zeros (n,n);
         [L,U,P,Delta] = LUfac (A);
         [x,Delta] = LUsub (L,U,P,y);
         while (t > tol) & (k < m)
            y = x/norm(x,inf);
            [x,Delta] = LUsub (L,U,P,y);
            c = dot(y,x)/dot(y,y);
            r = c*y - x;
            t = norm (r,inf);
            k = k + 1;
            if display
               fprintf ('\n%g & %.7f',k,c);
               fprintf (' & %.7f',y);
               fprintf (' & %.7f \\\\',t);
            end
         end  
         c = 1/c;
      end
   end   
      
% Finalize 

   x = y;
   if display
      fprintf ('\n');  
      wait;
   end
%----------------------------------------------------------------