function [U, S, V]=bidiagonalize(a)
% Naiv programmierte Bidiagonalisierung zur Vorbereitung einer SVD
% Am Ende: S bidiagonale obere Dreiecksmatrix,
% und U'*S*V=a mit U, V orthogonal
[m,n]=size(a);
U=eye(m);
V=eye(n);
mn=min(m-1,n-1);
for i=1:mn
    um=zeros(m,1);
    um(i:end)=householder(a(i:end,i));
    a=a-2*um*(um'*a);
    U=U-2*um*(um'*U);
    %
    un=zeros(n,1);
    un(i+1:end,1)=householder(a(i,i+1:end)');
    a=a-2*(a*un)*un';
    V=V-2*un*(un'*V); 
end
if m<n
    un=zeros(n,1);
    un(m+1:end,1)=householder(a(m,m+1:end)');
    a=a-2*(a*un)*un';
    V=V-2*un*(un'*V); 
end
if n<m
    um=zeros(m,1);
    um(n:end)=householder(a(n:end,n));
    a=a-2*um*(um'*a);
    U=U-2*um*(um'*U);
end  
S=a;
