Mercurial > matrix-functions
comparison toolbox/bandred.m @ 2:c124219d7bfa draft
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Thu, 07 May 2015 18:36:24 +0200 |
parents | 8f23314345f4 |
children |
comparison
equal
deleted
inserted
replaced
1:e471a92d17be | 2:c124219d7bfa |
---|---|
1 function A = bandred(A, kl, ku) | |
2 %BANDRED Band reduction by two-sided unitary transformations. | |
3 % B = BANDRED(A, KL, KU) is a matrix unitarily equivalent to A | |
4 % with lower bandwidth KL and upper bandwidth KU | |
5 % (i.e. B(i,j) = 0 if i > j+KL or j > i+KU). | |
6 % The reduction is performed using Householder transformations. | |
7 % If KU is omitted it defaults to KL. | |
8 | |
9 % Called by RANDSVD. | |
10 % This is a `standard' reduction. Cf. reduction to bidiagonal form | |
11 % prior to computing the SVD. This code is a little wasteful in that | |
12 % it computes certain elements which are immediately set to zero! | |
13 % | |
14 % Reference: | |
15 % G.H. Golub and C.F. Van Loan, Matrix Computations, second edition, | |
16 % Johns Hopkins University Press, Baltimore, Maryland, 1989. | |
17 % Section 5.4.3. | |
18 | |
19 if nargin == 2, ku = kl; end | |
20 | |
21 if kl == 0 & ku == 0 | |
22 error('You''ve asked for a diagonal matrix. In that case use the SVD!') | |
23 end | |
24 | |
25 % Check for special case where order of left/right transformations matters. | |
26 % Easiest approach is to work on the transpose, flipping back at the end. | |
27 flip = 0; | |
28 if ku == 0 | |
29 A = A'; | |
30 temp = kl; kl = ku; ku = temp; flip = 1; | |
31 end | |
32 | |
33 [m, n] = size(A); | |
34 | |
35 for j=1 : min( min(m,n), max(m-kl-1,n-ku-1) ) | |
36 | |
37 if j+kl+1 <= m | |
38 [v, beta] = house(A(j+kl:m,j)); | |
39 temp = A(j+kl:m,j:n); | |
40 A(j+kl:m,j:n) = temp - beta*v*(v'*temp); | |
41 A(j+kl+1:m,j) = zeros(m-j-kl,1); | |
42 end | |
43 | |
44 if j+ku+1 <= n | |
45 [v, beta] = house(A(j,j+ku:n)'); | |
46 temp = A(j:m,j+ku:n); | |
47 A(j:m,j+ku:n) = temp - beta*(temp*v)*v'; | |
48 A(j,j+ku+1:n) = zeros(1,n-j-ku); | |
49 end | |
50 | |
51 end | |
52 | |
53 if flip | |
54 A = A'; | |
55 end |