function [A,d]=HouseholderQR(A);
% HOUSEHOLDERQR computes the QR-decomposition of a matrix
%   [A,d]=HouseholderQR(A) computes an implicit QR-decomposition A=QR
%   using Householder transformations. The output matrix A contains
%   the Householder vectors u and the upper triangle of R. The
%   diagonal of R is stored in the vector d.

[m,n]=size(A);
for j=1:n,
  s=norm(A(j:m,j));
  if s==0, error('rank(A)<n'), end
  if A(j,j)>=0, d(j)=-s; else d(j)=s; end
  fak=sqrt(s*(s+abs(A(j,j))));
  A(j,j)=A(j,j)-d(j);
  A(j:m,j)=A(j:m,j)/fak;
  if j<n,       % transformation of the rest of the matrix G:=G-u*(u'*G)
    A(j:m,j+1:n)=A(j:m,j+1:n)-A(j:m,j)*(A(j:m,j)'*A(j:m,j+1:n));
  end
end
