function [U, Q, L, S] = gqr(A, B, partial)
%GQR Generalized QR factorization.
% [U, Q, L, S] = GQR(A, B, partial) factorizes
% the m-by-n A and p-by-n B, where m >= n >= p, as
% A = U*L*Q^T, B = S*Q^T, with U and Q orthogonal
% and L = [0; L1], S = [S1 0], with L1 and S1 lower triangular.
% If a nonzero third argument is present then only a partial reduction
% of A is performed: the first p columns of A are not reduced to
% triangular form (which is sufficient for solving the LSE problem).
% Reference:
% N. J. Higham, Accuracy and Stability of Numerical Algorithms,
% Second edition, Society for Industrial and Applied Mathematics,
% Philadelphia, PA, 2002; sec. 20.9.
[m n] = size(A);
[p n1] = size(B);
if nargin < 3, partial = 0; end
if n ~= n1, error('A and B must have same number of columns!'), end
if partial
limit = p+1;
else
limit = 1;
end
[Q, S] = qr(B');
S = S';
U = eye(m);
A = A*Q;
% QL factorization of AQ.
for i = n:-1:limit
% Vector-reversal so that Householder eliminates leading
% rather than trailing elements.
temp = A(1:m-n+i,i); temp = temp(end:-1:1);
[v, beta] = gallery('house',temp);
v = v(end:-1:1);
temp = A(1:m-n+i,1:i);
A(1:m-n+i,1:i) = temp - beta*v*(v'*temp);
% Put zeros where they're supposed to be!
A(1:m-n+i-1,i) = zeros(m-n+i-1,1);
temp = U(:,1:m-n+i);
U(:,1:m-n+i) = temp - beta*temp*v*v';
end
L = A;