householder1.m

来自「Solve linear system and test it.」· M 代码 · 共 30 行

function [R,v] = householder(A)
[m,n] = size(A);
if m>n,
NumberOfReflections = n;
else
NumberOfReflections = m - 1;
end
R = A;
v = cell(NumberOfReflections,1);
for k = 1:NumberOfReflections,
x = R(k:m,k);
xnorm = norm(x);
if xnorm>0,
% Compute the normal vector of the reflector
v{k} = -x;
if x(1)<0||x(1)>0,
    v{k}(1) = v{k}(1) - sign(x(1))*xnorm;
else
    v{k}(1)=(-1)*xnorm;
end
v{k} = (sqrt(2)/norm(v{k}))*v{k};
% Update R
for j=k:n,
    R(k:m,j)=R(k:m,j)-(v{k}'*R(k:m,j))*v{k}
end
else
v{k} = zeros(m-k+1,1);
end
end
v