function [q,e,A]=Bidiagonalize(A)
% BIDIAGONALIZE bidiagonalizes a matrix with Householder reflections
%   [q,e,A]=Bidiagonalize(A) computes B=diag(q)+diag(e(2:n),1) such
%   that A=P B Q' using Householder reflexions. A is overwritten with
%   the Householder-vectors. 

[m,n]=size(A); 
for i=1:n
  s=norm(A(i:m,i));                     % transform A(i:m,i) to 
                                        % (q_i,0,...,0)
  if s==0, q(i)=0;
  else 
    if A(i,i)>0, q(i)=-s; else q(i)=s; end 
    fak=sqrt(s*(s+abs(A(i,i)))); 
    A(i,i)=A(i,i)-q(i); 
    A(i:m,i)=A(i:m,i)/fak;
    A(i:m,i+1:n)=A(i:m,i+1:n)-A(i:m,i)*(A(i:m,i)'*A(i:m,i+1:n));
  end
  if i<n,
    s=norm(A(i,i+1:n));                 % tranformation A(i,i+1:n) to 
                                        % (e_i,0...0)
    if s==0, e(i)=0; 
    else
      if A(i,i+1)>0, e(i)=-s; else e(i)=s; end 
      fak=sqrt(s*(s+abs(A(i,i+1)))); 
      A(i,i+1)=A(i,i+1)-e(i);
      A(i,i+1:n)=A(i,i+1:n)/fak; 
      A(i+1:m,i+1:n)=A(i+1:m,i+1:n) - ...
                    (A(i+1:m,i+1:n)*A(i,i+1:n)')*A(i,i+1:n);
    end 
  end 
end                                     % insert 0 element in e (see
e=[0 e];                                % notation of bidiagonal matrix)
