function A = bandred(A, kl, ku)
%BANDRED  Band reduction by two-sided unitary transformations.
%   B = BANDRED(A, KL, KU) is a matrix unitarily equivalent to A
%   with lower bandwidth KL and upper bandwidth KU
%   (i.e. B(i,j) = 0 if i > j+KL or j > i+KU).
%   The reduction is performed using Householder transformations.
%   If KU is omitted it defaults to KL.

%   Called by RANDSVD.
%   This is a `standard' reduction.  Cf. reduction to bidiagonal form
%   prior to computing the SVD.  This code is a little wasteful in that
%   it computes certain elements which are immediately set to zero!
%
%   Reference:
%   [1] G. H. Golub and C. F. Van Loan, Matrix Computations, third
%       edition, Johns Hopkins University Press, Baltimore, Maryland,
%       1996, Sec. 5.4.3.
%
%   Nicholas J. Higham, Dec 1999.
%   Copyright 1984-2004 The MathWorks, Inc. 
%   $Revision: 1.4.4.2 $  $Date: 2004/01/24 09:21:39 $

if nargin == 2, ku = kl; end

if kl == 0 && ku == 0
   error('MATLAB:bandred:UseSVDforDiagonal',...
         'You''ve asked for a diagonal matrix.  In that case use the SVD!')
end

% Check for special case where order of left/right transformations matters.
% Easiest approach is to work on the transpose, flipping back at the end.
flip = 0;
if ku == 0
   A = A';
   temp = kl; kl = ku; ku = temp; flip = 1;
end

[m, n] = size(A);

for j=1 : min( min(m,n), max(m-kl-1,n-ku-1) )

    if j+kl+1 <= m
       [v, beta] = house(A(j+kl:m,j));
       temp = A(j+kl:m,j:n);
       A(j+kl:m,j:n) = temp - beta*v*(v'*temp);
       A(j+kl+1:m,j) = zeros(m-j-kl,1);
    end

    if j+ku+1 <= n
       [v, beta] = house(A(j,j+ku:n)');
       temp = A(j:m,j+ku:n);
       A(j:m,j+ku:n) = temp - beta*(temp*v)*v';
       A(j,j+ku+1:n) = zeros(1,n-j-ku);
    end

end

if flip
   A = A';
end