% Householder QR factorization
% Transform an m by n matrix with m >= n
% to upper triangular form.
% Algorithm 3.1 on p.124 Heath
% Example 3.8 on p.125 Heath
% HouseholderQR.m
% K. Ming Leung, 01/14/03 

clear all; format long g;
A = [ 1  0  0
      0  1  0
      0  0  1
     -1  1  0
     -1  0  1
      0 -1  1];
b = [ 1237
      1941
      2417
       711
      1177
       475];
   
[Q,R]=qr(A);    % Matlab's qr factorization
c=Q'*b;         % transform rt-hand side using Matlab

[m,n] = size(A); 
for k = 1:n      
    v = zeros(m,1);
    alpha = -sign(A(k,k))*norm(A(k:m,k));
    v(k:m) = A(k:m,k); v(k) = v(k) - alpha;
    beta = (v'*v)/2;
    if (beta~=0)
        for j = k:n
            gamma = (v'*A(:,j))/beta;
            A(:,j) = A(:,j) - gamma*v;
        end
        gamma = v'*b/beta;
        b = b - gamma*v;
    end
end
% A should be the same as R
A,
R,
% b should be the same as c
b,
c,
x=backsubstitution(A(1:n,1:n),b,n),