holderk.m

来自「基于Matlab语言编写的」· M 代码 · 共 25 行

function [H,u,p]=holderk(x,k)
%求householder矩阵Hk
n=length(x);
if k<1|k>=n
    disp('k值错误');
    return;
end
M=max(abs(x));
if M==0
    disp('M=0');
    return;
else
    x=x/M;
end
s=norm(x(k:n));
if x(k)<0
    s=-s;
end
x(k)=x(k)+s;
p=s*x(k);
u=zeros(n,1);
u(k:n)=x(k:n);
H=eye(n,n)-p\u*u';
u=M*u;
p=M*M*p;