function [R,Q]=rqfac(A)
% function [R,Q]=rqfac(A)
% returns the RQ factorization of A

[m,n]=size(A);

if (m~=n)
	error('A must be square')
end

Q=eye(m);

for i=n:-1:1
	u=house(A(i,1:i));
	if (1)
		% Matrix formulation.
		H=eye(length(u))-2*u'*u;
		A(1:i,1:i)=A(1:i,1:i)*H;
		Q(:,1:i)=Q(:,1:i)*H;
	else
		% Minimal flops.
		for j=1:i
			A(j,1:i)=housemul(u,A(j,1:i));
		end
		for j=1:m
			Q(j,1:i)=housemul(u,Q(j,1:i));
		end
	end
end

Q=Q';
R=A;

function u=house(x)
% function u=house(x)
% returns the Householder vector such that all elements but the last in 
% (I-2*u*u')*x is zero.

n=norm(x);
u=x;
u(end)=u(end)+sign(u(end))*n;
u=u/norm(u);

function z=housemul(u,y)
% function z=housemul(u,y)
% Applies the Householder transformation (I-2*u*u') to y from the right
% without actually calculating (I-2*u*u')

z=y-2*(y*u')*u;

