% Script to test HouseholderQR.m
% for performing Householder QR factorization
% to transform an m by n matrix A with m >= n
% to upper triangular form.
% Result gives an m by m orthogonal matrix QQ
% & an m by n upper triangular matrix RR 
% such that QQ*RR = A.
% Algorithm 3.1 on p.124 Heath
% Example 3.12 & 13 on p.135 & 136 Heath
% HouseholderQR.m
% K. Ming Leung, 04/6/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];
     
% % 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 rank defficient
%   its R is therefore nearly singular:
%   R =    [-1.2848416244814        -0.927442711455564
%                          0      0.000130261090298666
%                          0                         0]
%   

[Q,R]=qr(A);        % Matlab's qr factorization
c=Q'*b;             % transform rt-hand side using Matlab

[m,n]=size(A);
[QQ,RR]=HouseholderQR(A);
% Compare with Matlab results:
errorRR=max(max(abs(RR-R))),
errorQQ=max(max(abs(QQ-Q))),
% Solve the least squares problem:
cc = QQ.'*b;
x=backsubstitution(RR(1:n,1:n),cc,n),