% 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),