house.m

matlab算法集 matlab算法集

function [B,Q] = house (A)
%----------------------------------------------------------------------- 
% Usage:       [B,Q] = house (A);
%
% Description: Convert the n by n matrix A to upper-Hessenberg form
%              using the Householder transformations.
%
% Inputs:      A = n by n matrix
%
% Outputs:     B = n by n upper-Hessenberg version of A.  B is zero
%                  below the subdiagaon.  If A is symmetric, then
%                  B is tridiagonal.
%              Q = n by n orthogonal transformation matrix used to
%                  find B using B = Q'AQ.
%
% Notes:       house is used to as a preprocessing step to speed
%              the convergence of iterative methods
%----------------------------------------------------------------------- 
            
% Initialize, check elements below subdiagonal 

   beta = 0;
   [m,n] = size(A); 
   u = zeros (n,1);
   P = zeros (n,n);
   B = A;
   Q = eye(n,n);
   if n < 3 
      return;
   end
   for k = 1 : n-2
      for j = k+2 : n
         beta = beta + abs(B(k+2,k));
      end
   end
   if beta < eps			
      return;   
   end

% Perform Householder transformations    

   for k = 1 : n-2

% Compute alpha  

      alpha = 0;
      for j = k+1 : n
         alpha = alpha + B(j,k)*B(j,k);
      end
      alpha = sign(B(k+1,k))*sqrt(alpha);
      
% Compute u       
      
      u = zeros(n,1); 
      for j = k+1 : n
         u(j) = B(j,k);
      end
      u(k+1) = u(k+1) + alpha;

% Compute P 
      
      P = eye (n,n);
      beta = 0.5*dot(u,u);
      for i = 1 : n
         for j = 1 : n
            P(i,j) = P(i,j) - u(i)*u(j)/beta;
         end
      end
   
% Update Q and B 
            
      Q = Q*P;
      B = P*B*P;

   end