jacobi.m

matlab算法集 matlab算法集

function [c,P,r] = jacobi (A,tol,m)
% ---------------------------------------------------------------------- 
% Usage:       [c,P,r] = jacobi (A,tol,m)
%
% Description: Find the eigenvalues and eigenvectors of a symmetric
%              n by n matrix A using Jacobi's iterative method.
%
% Inputs:      A   = n by n matrix
%              tol = error tolerance used to terminate search (tol >= 0)
%              m   = maximum number of iterations (m >= 1)
%
% Outputs:     c   = n by 1 vectors containing eigenvalues ordered
%                    by decreasing magnitude
%              P   = n by n matrix containing eigenvectors.  The kth
%                    eigenvector is in the kth column of P
%              r   = number of iterations. If r < m, then A was
%                    successfully converted to a diagonal matrix
%                    with the elements off the diagonal having
%                    magnitudes less than tol.
% ---------------------------------------------------------------------- 

% Initialize

   r = 0;
   n = size(A,1);
   m = args (m,1,m,3,'jacobi');
   if norm(A-A',1) > 100*eps
      fprintf ('\nThe matrix passed to jacobi must be symmetic.\n');
      return
   end
   P = zeros(n,n);
   c = zeros(n,1);
   if n == 1
      c(1) = A(1,1);
      P(1,1) = 1;
      return;
   end
   P = eye(n,n);
   [B,Q] = house (A);
      
% Perform iterations to diagonalize b         

   [p,q] = getmax(B);
   while (abs(B(p,q)) >= tol) & (r < m) 

% Compute phi 

      phi = 0.5*atan2(2*B(p,q),B(p,p)-B(q,q));
      cphi = cos (phi);
      sphi = sin (phi);
      
% Compute P and B = BP       
      
      for i = 1 : n
      	 a = P(i,p);
      	 b = P(i,q);
      	 P(i,p) =  cphi*a + sphi*b;
      	 P(i,q) = -sphi*a + cphi*b;
      	 a = B(i,p);
      	 b = B(i,q);
      	 B(i,p) =  cphi*a + sphi*b;
      	 B(i,q) = -sphi*a + cphi*b;
      end	 

% Compute B = P'BP 
      	 
      for j = 1 : n
      	a = B(p,j);
      	b = B(q,j);
         B(p,j) =  cphi*a + sphi*b;
         B(q,j) = -sphi*a + cphi*b;
      end
	 
% Find largest element */

      [p,q] = getmax (B);
      r = r + 1;

   end
   
% Sort eigenvalues and eigenvectors 

   c = diag(B);
   P = Q*P;
   for i = 1 : n
      for j = i : n
         if abs(c(j)) > abs(c(i))
            a = c(i);
            c(i) = c(j);
            c(j) = a;
            d = P(:,i);
            P(:,i) = P(:,j);
            P(:,j) = d;
         end   
      end
   end
% ----------------------------------------------------------------------