function [q, a, p, rk]=qrpivot(a)
% Householderverfahren zur QR-Zerlegung
% mit Pivotisierung und Rangentscheid. 
% Inputmatrix a darf nichtquadratisch sein.
% Ergebnis: a(:,p(:))=q*r  mit Permutationsvektor p
% und Rang in rk
[m n]=size(a);
p=1:n; % Permutation der Spalten
q=eye(m); % wir wollen auch q=Q' berechnen
qsum=sum(a.*a); % Spaltenquadratsummen als Zeilenvektor
mn=min(m,n);   % Schrittzahl des Verfahrens
for k=1:mn   % Zerlegungsschleife
    % Maximum der laufenden Spaltenquadratsummen
    [qsmax kmax]=max(qsum(1,:)); 
    % Diese werden upgedatet, und 
    % fuer die schon orthogonalisierten
    % Spalten sollten die Quadratsummen etwa Null sein.
    if kmax<k    % wenn das Max. vorne ist, 
        % ist irgendwas faul, nicht weiterrechnen
        rk=k-1;
        break;
    end
    if qsmax<eps % ebenso, wenn keine brauchbare 
                 % Quadratsumme mehr da ist
        rk=k-1;
        break;
    end
    % Volle Spaltenvertauschung
    help=a(:,kmax);
    a(:,kmax)=a(:,k);
    a(:,k)=help;
    % Permutation veraendern
    hp=p(kmax);
    p(kmax)=p(k);
    p(k)=hp;
    % Quadratsumme auch vertauschen
    hp=qsum(kmax);
    qsum(kmax)=qsum(k);
    qsum(k)=hp;
    % So, jetzt koennen wir normal householden
    u=zeros(m,1);
    u(k:end,1)=householder(a(k:end,k));
    % wir rechnen wieder mal primitiv
    a=a-2*u*(u'*a);
    q=q-2*u*(u'*q);
    % aber wir muessen noch die Quadratsummen updaten 
    qsum=qsum-a(k,:).*a(k,:);
    rk=k;
end
q=q';
