function  A = bidiag(A)
% bidiag  Reduction to bidiagonal form.
%
%    B = bidiag(A) reduces the matrix A to bidiagonal form
%    using the function hzero2 to zero-out entries under the
%    diagonal and elements A(i,i+2), A(i,i+3), ..., A(i,n)

[m,n] = size(A);

n = min(m-1, n);

for i = 1:n
   A = hzero2(A,i,i);
   if i <= n-1
      A = hzero2(A,i,i+1,1);
   end
end
A = bid(diag(A),diag(A,1));
