function [x r2 Q R]=Qrdes_3(A,b)
% Resolución de Ax=b mediante transformaciones de Householder; calcula Q y R
%
 [m n]=size(A); x=zeros(n,1); Q=eye(m); 
 for j=1:n
    w=Housv(A(j:m,j));                             % Householder de a(j:m,j)
    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));
    Qk=eye(m); 
    Qk(j:m,j:m)=eye(m+1-j)-2*(w*w'); 
    Q=Qk*Q;
 end
 for i=n:-1:1                                      % Rx=b
    x(i)=(b(i)-x(i+1:n)'*A(i,i+1:n)')/A(i,i);
 end
 r2=norm(b(n+1:m))^2;                              % Residuos al cuadrado
 R=triu(A); Q=Q';                                  % Matrices R y Q
end

function w = Housv(x)
% Transformación de Householder del vector x.
 m=max(abs(x)); w=x/m;
 sw=1; if w(1)<0, sw=-1; end
 w(1)=w(1)+sw*norm(w);
 w=w/norm(w); w=w(:);
end
