--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/bandred.m	Thu May 07 18:36:24 2015 +0200
@@ -0,0 +1,55 @@
+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:
+%         G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
+%         Johns Hopkins University Press, Baltimore, Maryland, 1989.
+%         Section 5.4.3.
+
+if nargin == 2, ku = kl; end
+
+if kl == 0 & ku == 0
+   error('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