  function ShowGivensQR()
% Illustrate Givens QR factorization.
% Given mxn A, compute an mxm orthogonal matrix Q (as a product of Givens 
% rotations) so that A = Q*R where R is mxn and upper triangular

format short
m = 6;
n = 4;
A = randn(m,n);
Q = eye(m,m);
pause(2)
for j=1:n
    for i=m-1:-1:j
        % Determine c and s so [c s;-s c]'*A(i:i+1,j) = [*;0]
        [c,s] = Givens(A(i,j),A(i+1,j));
        clc
        % Update A and Q...
        A(i:i+1,j:n) = [c s;-s c]'*A(i:i+1,j:n)
        Q(:,i:i+1) = Q(:,i:i+1)*[c s;-s c]
        pause(2)
    end
end


function [c,s] = Givens(a,b)
% function [c,s] = Givens(a,b)
% Givens rotation computation
% Determines cosine-sine pair (c,s) so that [c s;-s c]'*[a;b] = [r;0]
if b==0
    c = 1; s = 0;
else
    if abs(b)>abs(a)
        tau = -a/b; s = 1/sqrt(1+tau^2); c = s*tau;
    else
        tau = -b/a; c = 1/sqrt(1+tau^2); s = c*tau;
    end
end
% This is more reliable than setting c = a/sqrt(a^2+b^2) and 
% s = -b/sqrt(a^2+b^2). 