% Householder QR factorization
% Transform an n by n SQUARE matrix A
% to upper triangular form.
% Algorithm 3.1 on p.124 Heath
% Example 3.8 on p.125 Heath
% HouseholderQRsquare.m
% K. Ming Leung, 01/14/03
% Example 7.9: Fausett, p.204. 

clear all;
A = [21  7  -1
      5  7   7
      4 -4  20]; 
% Resulting A is the upper triangular matrix R:
% A=[-21.9545   -7.5611   -4.2816
%           0   -7.5386    9.3351
%           0    0.0000   18.5613];


[Q,R]=qr(A),
[n,n] = size(A); v = zeros(n,1);
for k = 1:n      
    alpha = -sign(A(k,k))*norm(A(k:n,k));
    v(k+1:n) = A(k+1:n,k);      v(k) = A(k,k) - alpha;
    beta = (v(k:n)'*v(k:n))/2;
    if (beta~=0)
        gamma = (v(k:n)'*A(k:n,k))/beta;
        A(k,k) = A(k,k) - gamma*v(k);
        A(k+1:n,k) = 0;
        for j = k+1:n
            gamma = (v(k:n)'*A(k:n,j))/beta;
            A(k:n,j) = A(k:n,j) - gamma*v(k:n);
        end
    end
end
RR=A, 
maxError = max(max(abs(R-RR))),