Mercurial > octave-nkf
changeset 20331:0b9d23557506
gallery: fix randsvd by adding missing dependency qmult().
* scripts/special-matrix/gallery.m (randsvd) was copied from the
Test Matrix toolbox by Nicholas J. Higham. It made use of qmult()
which was also part of that toolbox but was left behind. This qmult()
implementation is also recovered from the Test Matrix toolbox. Note
that Octave itself used to have qmult() which is part of the legacy
quaternion package (now also removed from the new quaternion package).
See cset 21904fe299c8 for when qmult() was removed from Octave.
Also fix the default value for KL and KU so that a 2 element vector
can be used as N.
| author | Carnë Draug <carandraug@octave.org> |
|---|---|
| date | Fri, 03 Jul 2015 16:18:33 +0100 |
| parents | 557979395ca9 |
| children | 26fc9bbb8762 |
| files | scripts/special-matrix/gallery.m |
| diffstat | 1 files changed, 49 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/special-matrix/gallery.m Fri Jul 03 15:01:38 2015 +0100 +++ b/scripts/special-matrix/gallery.m Fri Jul 03 16:18:33 2015 +0100 @@ -2150,7 +2150,8 @@ endfunction -function A = randsvd (n, kappa = sqrt (1/eps), mode = 3, kl = n-1, ku = kl) +function A = randsvd (n, kappa = sqrt (1/eps), mode = 3, kl = max (n) -1, + ku = kl) ## RANDSVD Random matrix with pre-assigned singular values. ## RANDSVD(N, KAPPA, MODE, KL, KU) is a (banded) random matrix of order N ## with COND(A) = KAPPA and singular values from the distribution MODE. @@ -2851,6 +2852,53 @@ endif endfunction +## NOTE: qmult is part of the Test Matrix Toolbox and is used by randsvd() +function B = qmult (A) + ## QMULT Pre-multiply by random orthogonal matrix. + ## QMULT(A) is Q*A where Q is a random real orthogonal matrix from + ## the Haar distribution, of dimension the number of rows in A. + ## Special case: if A is a scalar then QMULT(A) is the same as + ## QMULT(EYE(A)). + ## + ## Called by RANDSVD. + ## + ## Reference: + ## G.W. Stewart, The efficient generation of random + ## orthogonal matrices with an application to condition estimators, + ## SIAM J. Numer. Anal., 17 (1980), 403-409. + + [n, m] = size (A); + + ## Handle scalar A + if (isscalar (A)) + n = A; + A = eye (n); + endif + + d = zeros (n); + + for k = n-1:-1:1 + ## Generate random Householder transformation. + x = randn (n-k+1, 1); + s = norm (x); + sgn = sign (x(1)) + (x(1) == 0); # Modification for sign(1)=1. + s = sgn*s; + d(k) = -sgn; + x(1) = x(1) + s; + beta = s*x(1); + + ## Apply the transformation to A. + y = x'*A(k:n,:); + A(k:n,:) = A(k:n,:) - x*(y/beta); + endfor + + ## Tidy up signs + for i = 1:n-1 + A(i,:) = d(i)*A(i,:); + endfor + A(n,:) = A(n,:) * sign (randn); + B = A; +endfunction ## BIST testing for just a few functions to verify that the main gallery ## dispatch function works.