gusucode.com > elmat工具箱matlab源码程序 > elmat/private/bandred.m

    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, third
%   edition, Johns Hopkins University Press, Baltimore, Maryland,
%   1996, Sec. 5.4.3.
%
%   Nicholas J. Higham
%   Copyright 1984-2013 The MathWorks, Inc.

if nargin == 2
    ku = kl;
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);
niter = min( min(m,n), max(m-kl-1,n-ku-1) );

for j = 1:niter
    
    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) = 0;
    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) = 0;
    end
    
end

if flip
    A = A';
end