function [r,q] = qr_house(x)

% purpose: 用Householder鏡像變換計算QR分解，x = qr，size(x)=[m,n]
% 目前只適用于 m ≥ n 的情況，且size(q) = size(x)
% Data: 31/3/2008
% Version 1.0
% Copyright (C) 2008, Yang Chun
% 
% 輸入 - x 為列滿秩矩陣
% 輸出 - q 為x的正交基
%      - r 為正交基的坐標，為上三角矩陣
% 
% r = qr_house(x) 返回正交基坐標
% [r,q] = qr_house(x) 返回正交基與正交基坐標

[m,n] = size(x);
if m > n
    iterend = n;
else
    iterend = n - 1;
end
for k1 = 1:iterend
    [v,beta] = house(x(k1:m,k1));
    G = v'*x(k1:m,k1:n);
    S = beta*v;
    x(k1:m,k1:n) = x(k1:m,k1:n) - S*G;
    x(k1+1:m,k1) = v(2:m-k1+1);
    d(k1) = beta;
end
clear G S

if nargout == 1 
    r = triu(x(1:n,:));
else
    r = triu(x(1:n,:));
    q = eye(m,n);
    for k1 = iterend:-1:1
        v = [1;x(k1+1:m,k1)];
        G = v'*q(k1:m,k1:n);
        S = d(k1)*v;
        q(k1:m,k1:n) = q(k1:m,k1:n) - S*G;
    end
end