Mercurial > matrix-functions
comparison toolbox/bandred.m @ 0:8f23314345f4 draft
Create local repository for matrix toolboxes. Step #0 done.
| author | Antonio Pino Robles <data.script93@gmail.com> |
|---|---|
| date | Wed, 06 May 2015 14:56:53 +0200 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:8f23314345f4 |
|---|---|
| 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 |