function givensqr_demo(A)
%GIVENSQR QR factorization by Givens rotation 
%[Q,R] = givqr(A) returns an orthogonal matrix Q
%and a matrix R of the same size as A 
%with zeros below the diagonal, such that A = QR
%Input: m x n matrix A
%Output: m x m orthogonal matrix Q and m x n
% upper triangular matrix R such that A = QR

[m,n] = size(A);

%Q will contain the product, in reverse order, of the
%inverses (transposes) of the orthogonal Givens matrices
%that transform A into an upper triangular matrix
Q = eye(m,m);

disp('Initial matrix:');
A
spy(A);
pause;

for i= 1:min(m-1,n)
	fprintf('Beginning to zero-out all elements below (%d,%d)\n', i,i);
	for j = i+1:m
		%get the Givens parameters
		[c,s] = givensparms(A(i,i),A(j,i));
		fprintf('Rotation J(%d,%d,%f,%f) is\n',i,j,c,s);
		J = eye(m);
		J(i,i) = c;
		J(j,j) = c;
		J(j,i) = -s;
		J(i,j) = s;
		J
		disp('Apply the rotation to A.');
		%get the current Givens matrix
		%apply the rotation to A
		A = givensmul(A,i,j,c,s);
		A
		%maintain the product of the Givens matrices
		Q = givensmul(Q,i,j,c,s);
		disp('Current product of Givens rotations');
		Q

		% set almost-zero entries to true zero
		% display matrix and wait for a key
% 		A(find(abs(A)<1.e-13)) = 0;
		A(abs(A)<1.e-13) = 0;
		spy(A);
		pause;

		disp('======================================');
	end
	disp('____________________________________________________________');
end 
R = A
Q = Q'
