function [Q R] = qrdecomp(A)
% function [Q R] = qrdecomp(A)
% Bestimmt die QR-Zerlegung der Matrix A unter Verwendung von
% Householdertransformationen
%
% Input:
% A...rechteckige Matrix
%
% Output:
% Q...orthogonale Matrix
% R...rechte obere Dreiecksmatrix
%
% Es gilt:
% Q*R=A
%
%
% Wien, 28.04.2011
% Markus Grasmair



[m,n] = size(A);

% Initialisierung
R = A;
Q = eye(m);

for i=1:min(n,m)
    % Bestimmung des Vektors, entlang dem gespiegelt wird
    if(norm(R(i:end,i))==0)
        error('Die Matrix A hat nicht vollen Rang');
    end
    v = R(i:end,i)/norm(R(i:end,i));
    % der Befehl sign kann nicht verwendet werden, weil sign(0) = 0, wir
    % aber auch fuer diesen Fall als Ergebnis entweder +1 oder -1 brauchen
    if(v(1)>=0)
        v(1) = v(1) + 1;
    else
        v(1) = v(1) - 1;
    end
    
    % Anwendung der Householdertransformation auf R
    beta = 2/norm(v)^2;
    w = beta*(R(i:end,i:end)')*v;
    R(i:end,i:end) = R(i:end,i:end)-v*w';
    
    % Anwendung der Householdertransformation auf Q'
    wq = beta*(Q(i:end,1:end)')*v;
    Q(i:end,1:end) = Q(i:end,1:end)-v*wq';
end

Q = Q';