Implementazione in Matlab per risolvere $QRx=b$

Come sempre ci riportiamo al nostro problema di partenza, dunque alla risoluzione del sistema lineare $Ax=b$ tramite fattorizzazione $QR$. Dal momento che $Q^TA=R$ possiamo scrivere $Q^TAx=Q^Tb$ e dunque $Rx=\hat{b}$; il vettore $\hat{b}$ è il vettore dei termini noti modificato applicando la matrice $Q^T$, ma questa modifica può essere fatta mentre si esegue l'algoritmo, moltiplicando anch'esso per le matrici di Householder:

function x=solveQR(A,b)
%SOLVEQR Risolve il sistema lineare Ax=b fattorizzando la matrice
%   A come QR ed infine risolvendo il sistema
%                Rx=b^
%   dove b^ è il vettore dei termini noti aggiornato come b^=Q^Tb
%
%   x=SOLVEQR(A,b)
%
%   I parametri della funzione sono:
%       A -> la matrice dei coefficienti del sistema lineare
%       b -> il vettore dei termini noti
%
%   I valori di ritorno sono:
%       x -> il vettore soluzione del sistema lineare
%
%   See Also FATTQR
  [m,n]=size(A);
  x=b;
  for i=1:n
     alpha=norm(A(i:m,i));
     if alpha==0
        disp('No rango MAX!!');
        return
     end
     v1=A(i,i);
     if v1>=0
        v1=v1+alpha;
        s=1;
        A(i,i)=-alpha;
     else
        v1=v1-alpha;
        s=-1;
        A(i,i)=alpha;
     end
     A(i+1:m,i)=A(i+1:m,i)/v1;
     vt=[1;A(i+1:m,i)];
     beta=s*v1/alpha;
     A(i:m,i+1:n)=A(i:m,i+1:n)-(beta*[1;A(i+1:m,i)])*
                      ([1 A(i+1:m,i)']*A(i:m,i+1:n));
     b(i:n)=b(i:n)-(beta*[1 A(i+1:m,i)']*b(i:n))*[1;A(i+1:m,i)];
  end
  x=solveUT(A,b);