% Householder QR factorization
% Transform an m by n matrix A with m >= n
% to upper triangular form.
% Algorithm 3.1 on p.124 Heath
% Example 3.12 & 13 on p.135 & 136 Heath
% HouseholderQR2.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];
  
% Example from p.282 Kincaid & Cheney
A = [  63   41  -88
       42   60   51
        0  -28   56
      126   82  -71];
   
%   A = [ -1  1 0
%         -1  0 1
%          0 -1 1];   % Example 3.12
% This A above is rank deficient
% its R is therefore singular:

%   A=[ 0.913 0.659
%       0.780 0.563
%       0.457 0.330];   % Example 3.13
%   This above A is nearly randk defficient
%   its R is therefore nearly singular:
%   R =    [-1.2848416244814        -0.927442711455564
%                          0      0.000130261090298666
%                          0                         0]
%   
% 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); QQ=eye(m);
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;
    H=eye(m)-v*v.'/beta;
    QQ = QQ*H;
    if (beta~=0)
        % for j=k
        gamma = (v'*A(:,k))/beta;
        A(k,k) = A(k,k) - gamma*v(k);
        A(k+1:m,k) = zeros(m-k,1);
        % for the rest of the columns
        for j = k+1: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
RR=A,
R,
QQ,
Q,
% b should be the same as c
% b,
% c,
% x=backsubstitution(A(1:n,1:n),b,n),