hessenberg约化matlab代码

以下是Hessenberg约化的MATLAB代码实现：

function [H,Q] = hessenberg(A)
% Hessenberg reduction of a square matrix A
% Input: A (n x n) - nonsingular matrix
% Output: H (n x n) - upper Hessenberg matrix
%         Q (n x n) - orthogonal matrix such that Q'*A*Q=H

n = size(A,1);
Q = eye(n);

for k = 1:(n-2)
    % Zero subdiagonal entries
    x = A((k+1):n,k);
    v = [sign(x(1))*norm(x);x];
    v(1) = v(1) + x(1);
    v = v/norm(v);
    A((k+1):n,k:n) = A((k+1):n,k:n) - 2*v*(v'*A((k+1):n,k:n));
    A(1:n,(k+1):n) = A(1:n,(k+1):n) - 2*(A(1:n,(k+1):n)*v)*v';
    Q(:,(k+1):n) = Q(:,(k+1):n) - 2*(Q(:,(k+1):n)*v)*v';
end

%Hessenberg matrix
H = triu(A, -1);
end

该代码输入为一个$n\times n$的非奇异矩阵$A$，输出为一个$n\times n$的上Hessenberg矩阵$H$以及一个$n\times n$的正交矩阵$Q$，满足$Q^T A Q = H$。