lin_solve.m

Solve linear system and test it.

function x = lin_solve(A,b)
[m,n] = size(A);
[R,v] = householder(A);
y = Qt_times_b(v,b,m);
x = R\y;
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)=-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

function y = Qt_times_b(v,b,m)
NumberOfReflections = length(v);
y = b;
for k =1:NumberOfReflections,
% Update y by multiplying Qk to y
y(k:m)=y(k:m)-(v{k}'*y(k:m))*v{k};
end