e351.m

来自「matlab算法集 matlab算法集」· M 代码 · 共 42 行

%-----------------------------------------------------------------------
% Example 3.5.1: Householder Transformation
%-----------------------------------------------------------------------
   clc                              
   clear                            
   n = 3;
   k = 1;
   alpha = 0;
   u = zeros (n,1);
   A = [1 -2 3; -2 4 1; 3 1 2];
   I = eye(n,n);

   fprintf ('Example 3.5.1: Househoulder Transformation\n');
   [B,P] = house (A);
   show ('A',A)
   show ('B',B)
   for j = k+1 : n
      alpha = alpha + A(j,k)*A(j,k);
   end
   alpha = sign(A(k+1,k))*sqrt(alpha);
   show ('alpha',alpha);
   for j = 1 : n
      if j < (k+1)
         u(j) = 0;
      elseif j == (k+1)
         u(j) = A(j,k) + alpha;
      elseif j >  (k+1)
         u(j) = A(j,k);
      end
   end
   show ('u',u)
   beta = dot(u,u)/2;
   for i = 1 : n
      for j = 1 : n
         P(i,j) = I(i,j) - u(i)*u(j)/beta;
      end
   end
   show ('P',P)
   show ('P''AP',P'*A*P)