% Taken directly from Ray Spiteri's notes,
% http://www.cs.usask.ca/~spiteri/M313/notes/Lecture16.pdf
% which themselves were taken from Trefethen and Bau:
R=triu(randn(50)+1i*randn(50)); % Set R to a 50x50 upper-triangular matrix
% with normal random entries
[Q,X]=qr(randn(50)+1i*randn(50)); % Set Q to a 50x50 random orthogonal matrix
% by orthogonalizing a random matrix
A=Q*R; % Set A=QR, up to rounding errors
[Q2,R2]=qr(A); % Compute QR factorization by Householder
%-------------------------------------
% Modified from here (RMC)
% Make sure the signs are the same even
% in the complex case (which this isn't now, but might be)
D1 = diag(diag(sign(R)));
D2 = diag(diag(sign(R2)));
% in each case Dbar * D = eye
% A = Q2*R2
% Good factoring initially
norm( A - Q2*R2, inf )
% A = Q2*Dbar*D*R2 
Dbar = conj( D1*D2 );
D    = D1*D2;
Q3 = Q2*Dbar;
R3 = D*R2;
norm( A - Q3*R3, inf )