Matlab - Error householder

function [x r res]=Mincua_QR(A,b)

% Resolución de Ax=b general mediante transformaciones de Householder

% Posible rango incompleto r

[m n]=size(A); x=zeros(n,1); tol=sqrt(eps); W=zeros(n,m); ip=1:n; r=n;

for j=1:n

jm=j; c=0;

for k=j:n

h=norm(A(j:m,k));

if h>c, c=h; jm=k; end

end

if jm~=j, A(:,[j jm])=A(:,[jm j]); ip([j jm]) = ip([jm j]); end

if j==m, break, end

w=Housv(A(j:m,j)); % Householder de A(j:m,j); luego a A y b

A(j:m,j:n)=A(j:m,j:n)-2*w*(w'*A(j:m,j:n)); b(j:m)=b(j:m)-2*w*(w'*b(j:m));

end

for j=1:n % Ver rango

if abs(A(j,j))<=tol, r=j-1; break, end

end

res=norm(b(r+1:m))^2; w1=zeros(r,n-r+1);

W(1:n,1:r)=A(1:r,1:n)'; % Trasp. de A = W

if r~=n

for i=r:-1:1

w1(i,1:n-r+1)=Housv([W(i,i);W(r+1:n,i)]); % Householder hacia W

W([i r+1:n],i:-1:1)=W([i r+1:n],i:-1:1)-2*w1(i,:)'*(w1(i,:)*W([i r+1:n],i:-1:1));

end

end

for i=r:-1:1 % Resol. Wx=c

x(i)=(b(i)-x(i+1:r)'*W(i+1:r,i))/W(i,i);

end

if r~=n

x(r+1:n)=0; % Aplicar a x ultimas Householder hacia W

for i=1:r, x([i r+1:n])=x([i r+1:n])-2*w1(i,:)'*(w1(i,:)*x([i r+1:n])); end

end

x(ip)=x(:); % Deshacer perm. col.

>> [x r res]=Mincua_QR

A =

     1     1     1

     1     2     4

     1     3     9

     1     4    16

b =

     2

     3

     5

     6

x =

    0.5000

    1.4000

    0.0000

r =

     3

res =

    0.2000